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Abstract 

The  goal  of  this  research  was  to  develop  a  realizable 
proportional-plus-integral  (PI)  feedback  tracker  to  control 
a  neutral  particle  beam.  The  design  is  based  on  detecting 
the  photo-electron  events  that  are  emitted  from  a  laser- 
excited  particle  beam  and  the  observed  events  are  used  by  a 
Meer  filter  to  locate  the  beam's  centerline.  The  observed 
events  are  modeled  by  a  Poisson  space-time  process  and  are 
composed  of  both  signal-  and  noise-induced  events.  The  Meer 
filter  is  a  stochastic  multiple  model  adaptive  estimator 
which  is  composed  of  a  bank  of  Snyder-Fishman  filters  and  is 
designed  to  distinguish  the  signal-induced  events  from  the 
noise-induced  events.  A  target  model  is  developed  from  a 
Gauss-Markov  acceleration  process,  and  the  target  states  are 
estimated  by  a  Kalman  filter.  The  "optimal"  PI  controller 
design  is  based  on  the  linear  quadratic  (LQ)  controller 
synthesis  technique  and  the  "assumed”  certainty  equivalence 
property,  and  the  Kalman  filter  provides  the  reference 
(target)  states  while  the  Meer  filter  supplies  controlled 
(beam)  states.  The  objectives  of  the  research  were  to  (1) 
select  the  "best”  cost  weighting  matrices  that  minimize  the 
RMS  tracker  error  and  enhance  robustness,  (2)  simplify  the 
Meer  filter  for  easier  on-line  usage,  (3)  complete  full- 

xiii 


scale  sensitivity  and  robustness  analyses  over  all  the 
Kalman  and  Meer  filter  parameters,  and  (4)  develop  on-line 
adaptive  estimation  of  those  parameters  that  greatly  affect 
stability  robustness  and  tracker  performance.  During  the 
research,  an  apparent  stability  problem  was  uncovered,  and 
fifth  objective  was  to  identify  the  source  of  the 
instability,  and  to  propose  a  solution  that  would  insure 
stability  during  parameter  variations. 
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This  research  is  motivated  in  part  by  the  problems  in 
neutral  particle  beam  pointing  and  tracking  currently  being 
investigated  at  the  Air  Force  Weapons  Laboratory,  and  at  the 
Rome  Air  Development  Center.  Their  goal  is  not  only  to 
estimate  the  position  and  direction  of  the  beam,  but  use 
that  information  in  an  optimal  way  to  control  the  pointing 
of  the  beam. 

A  method  for  sensing  the  location  has  been  proposed  in 
which  the  beam  is  illuminated  by  one  or  more  lasers  [23] . 

At  certain  angles  of  intersection  and  at  different  particle 
velocities,  the  particle  electrons  absorb  photons  from  the 
laser  beam  and  attain  a  higher  energy  state.  When  the 
particle  beam  electrons  relax  and  spontaneously  decay  to 
their  ground  energy  state,  they  expend  the  energy  as  light. 
By  detecting  the  light  energy,  the  position  of  the  beam  can 
be  inferred.  If  the  light  energy  were  to  arrive  at  the 
photo-detector  at  a  sufficient  signal  rate  as  to  produce  an 
observable  current,  it  might  be  modeled  by  a  continuous¬ 
time,  Gaussian  process  as  done  in  many  communication  and 
control  type  problems.  But,  the  assumption  upon  which  this 
thesis  is  based  is  that  the  photon  events  do  not  arrive  at 
such  a  sufficient  rate.  Instead,  a  discrete,  Poisson  space- 
time  process  is  used  to  describe  the  arrival  of  the 
individual  signal-induced  events  (the  photons)  and  the 
noise-induced  events  (caused  by  dark  currents  within  the 


detector,  or  other  outside  sources  of  noise) .  A  space-time 
point  process  is  a  stochastic  process  having  as  realization 
points  with  random  coordinates  in  both  time  and  space.  [9] . 
For  this  application,  the  time  between  events  will  assume  a 
conditional  Poisson  process  composed  of  Gaussian  spatially 
distributed  signal-induced  events,  and  uniformly  distributed 
noise-induced  events. 

1.1 

In  1975,  Snyder  and  Fishman  [25]  developed  an 
estimation  algorithm  called  the  Snyder-Fishman  filter.  The 
filter  is  a  minimum-mean-square  estimator  that  for  this 
application,  estimates  the  position  of  maximum  intensity 
from  the  arrival  of  the  signal-induced  events  (all  events 
are  assumed  to  be  signal  induced) .  The  filter  equations 
appear  very  similar  to  the  Kalman  filter  equations  with  the 
notable  exception  that  the  Snyder-Fishman  filter  is  based  on 
the  Poisson  space-time  process  and  thus  the  arrival  times  of 
the  events  are  not  known  a  priori.  Instead,  the  events 
arrive  as  part  of  the  Poisson  distributed  process.  In  the 
absence  of  noise-induced  events,  the  Snyder-Fishman  filter 
provides  good  results,  but  Santiago  [24]  found  the  filter's 
performance  is  severely  degraded  in  the  presence  of  noise. 

In  1982,  Meer  [18]  developed  an  adaptive  filter 
designed  to  estimate  the  position  of  maximum  intensity  from 
the  arrival  of  signal-induced  events  that  are  corrupted  by  a 
statistically  independent  noise  process.  This  was 
accomplished  through  a  multiple  model  structure  in  which  a 


bank  of  elemental  Snyder -Fishman  filters  are  based  on 
different  hypothesis  sequences.  The  hypotheses  define  which 
observed  events  were  assumed  to  be  due  to  the  signal  process 
and  which  were  assumed  to  be  due  to  the  noise  process,  and 
each  elemental  filter  is  allowed  to  process  the  observed 
events  only  when  its  associated  hypothesis  defined  the  event 
as  being  signal-induced.  Meer  was  the  first  to  apply  the 
point  process  to  the  neutral  particle  beam  problem,  and  the 
"Meer"  filter  outperformed  the  Snyder-Fishman  filter  in  a 
noisy  environment.  This  successful  filter  development  was 
the  primary  breakthrough  required  in  controlling  the  beam 
location. 

In  1983,  Zicker  [27]  conducted  a  feasibility  study  of  a 
simple  proportional  gain  controller,  considering  both  a 
regulator  design  to  null  out  variations  in  the  beam's 
location  and  a  tracker  design  to  maintain  the  beam  on  a 
maneuvering  target.  Zicker  synthesized  his  stochastic 
controller  designs  from  a  deterministic  optimal  LQ 
controller  assuming  full  state  feedback.  An  LQ  controller 
is  a  controller  design  based  on  a  linear  system  model  which 
uses  a  quadratic  performance  index  to  define  optimal 
control.  Then,  the  states  were  replaced  by  their  best 
estimates  according  to  the  principle  of  assumed  certain 
equivalence  [14] .  The  Meer  filter  was  used  to  provide  the 
beam  state  estimate  while  a  Kalman  filter  was  used  to 
provided  the  target  state  estimate. 


Over  the  next  two  years.  Moose  [20]  and  Jamerson  [8] 


replaced  Zicker's  proportion  gain  controller  with 
proportional-plus-integral  (PI)  controllers,  and  the  effort 
began  to  move  from  a  question  of  feasibility  to 
realizability  and  performance  potential.  The  PI  controller 
was  selected  because  it  possessed  several  characteristics 
which  make  it  ideally  suited  for  the  tracking  problem. 

First,  control  is  based  on  both  the  current  as  well  as  the 
integral  of  previous  tracking  errors,  thus  providing  a  type- 
one  system.  Second,  type-one  systems  are  able  to  reject 
constant  unmodeled  disturbances  that  arise  from  linearized 
models,  and  are  able  to  handle  non-zero  setpoint  control 
with  zero  steady  state  error.  Although  performance  improved 
considerably,  especially  when  the  controller  was  confronted 
with  unmodeled  disturbances,  it  came  at  the  expense  of  some 
additional  sluggishness  due  to  the  additional  integrator  in 
the  controller.  Also,  Jamerson  replaced  Zicker's  simple 
target  model,  which  was  based  on  a  Gauss-Markov  position 
model,  with  a  more  realistic  Gauss-Markov  acceleration 
model . 

1*2  9felg£&jygg 

Much  of  the  previous  work  was  spent  on  developing  a 
stochastic  filter  based  on  a  Poisson  space-time  point 
process  which  models  both  discrete  signal-  and  noise-induced 
events,  and  then  demonstrating  that  it  can  be  used  for 


which  has  the  stability  robustness  to  handle  "real  world" 
plant  variations.  Therefore,  the  objectives  of  this  thesis 
are:  (1)  to  evaluate  the  quadratic  cost  function  in  the  LQ 

synthesis  of  Jamerson's  PI  controller  and  select  the  best 
set  of  cost  weighting  matrices  that  minimize  the  RMS  tracker 
error  and  enhances  robustness,  (2)  to  evaluate  the 
possibilities  of  a  simplified  filter  to  replace  the  multiple 
model  structure  of  the  Meer  filter  in  order  to  reduce  the 


computational  loading,  {3)  to  complete  a  full-scale 
sensitivity  and  robustness  analysis  on  the  six  beam 
parameters  and  the  three  target  parameters  within  the  Meer 
and  Kalman  filters,  and  (4)  to  develop  on-line  adaptive 
estimation  for  those  parameters  that  greatly  affect 
stability  robustness  and  tracker  performance.  During  the 
research,  an  apparent  stability  problem  has  been  uncovered, 
and  a  fifth  objective  has  become  to  identify  the  source  of 
the  instability  and  to  propose  a  solution  that  insured 
stability  during  parameter  variations.  The  initial 
objective  of  evaluating  an  alternate  cost  weighting 
technique,  such  as  implicit  model  following  was  dropped  when 
the  apparent  stability  problem  surfaced.  Discussions  that 
pertain  to  developing  an  alternate  cost  weighting  technique 
remain  within  the  thesis  (see  Appendix  A)  as  an  aid  for 
future  research. 


An  overview  of  the  tracker  design  that  is  used 
throughout  the  thesis  is  shown  in  Figure  1-1.  The  purpose 
of  a  tracker  is  to  generate  a  control  input,  u(ti ) ,  that 
minimizes  the  difference  between  the  the  controller 
variable,  y,c  (t)  ,  and  the  reference  (i.e.,  target)  variable, 
Yr (t) .  In  other  words,  we  want  the  particle  beam  to  track 
the  target  position. 

The  target  model  is  generated  through  a  shaping  filter 
which  uses  an  exponentially  time-correlated  noise  process  to 
simulate  the  target  acceleration.  The  target  velocity  and 
position  states  are  generated  by  integrating  the 
acceleration  state  with  respect  to  time.  The  target  sensor 
is  modeled  as  a  sample  and  hold  device  that  obtains  a  noise- 
corrupted  measurement  of  the  target's  position,  £t  (ti  ) ,  at  a 
regular,  prespecifed  sample  rate.  The  measurements  are  used 
by  the  Kalman  filter  to  generate  a  target  state  estimate, 

St  • 

The  particle  beam  dynamics  (also  referred  to  as  the 
plant  dynamics)  are  modeled  as  an  exponentially  time- 
correlated  position  process,  the  output  of  a  shaping  filter 
having  one  dominate  pole.  The  position  of  the  beam  is 
inferred  through  the  observations  detected  at  the  surface  of 
the  photo-detector.  These  observations,  2.b  (t)  ,  are  assumed 
to  be  well  modeled  by  a  Poisson  space-time  point  process, 
and  can  be  of  either  signal  or  noise  origin.  The  Meer 
filter  is  designed  to  discriminate  between  the  two  types  of 


Figure  1-1.  System  Overview 


events,  and  generate  a  beam  state  estimate,  &b  (ti  )  ,  from  the 
signal-induced  events. 

The  target  state  estimate  and  the  beam  state  estimate 
are  fed  into  the  discrete-time  controller  algorithm  to 
generate  a  command  input  to  the  plant.  Because  the  target 
and  the  beam  processes  are  assumed  to  be  well  modeled  as 
Gaussian  processes  as  produced  by  linear  shaping  filter, 
"near"  optimal  controller  designs  can  be  developed  using  the 
linear  quadratic  (LQ)  full-state  feedback  controller 
synthesis  process,  and  the  assumed  certainty  equivalence 
property  for  incorporating  filters  into  the  loop.  The 
proportional  gain  and  proportional-plus-integral  controllers 
are  developed  in  theory  (presented  in  Chapter  3) ,  and  the  PI 
controller  is  implemented  and  evaluated  on  its  own  merit  and 
within  an  adaptive  multiple  model  controller  structure. 

The  performance  of  the  tracker  is  evaluated  by 
calculating  the  tracker  error,  e(t) ,  which  is  defined  as  the 
difference  between  the  target  and  beam  positions.  Because 
of  the  adaptive  nature  of  the  controller  and  the  time- 
variant  nature  of  the  Poisson  space-time  point  process,  the 
tracker  error  statistics  have  to  be  generated  with  a  Monte 
Carlo  simulation,  as  opposed  to  a  covariance  analysis  or 
similar  analytical  methods. 

1 . 4  Approach 

The  filter  and  controller  theory  is  developed  for  the 
n-dimensional  space,  but  for  simplicity,  the  specific 
tracking  system  of  interest  will  be  designed  and  evaluated 


in  one-space.  This  reduction  in  physical  dimensionality 
shall  simplify  the  physical  insights  about  the  design 
without  loss  of  generality.  Chapter  2  describes  the  Meer 
filter  in  considerable  detail.  The  tracker  model  and 
controller  designs  are  developed  in  Chapter  3.  This 
includes  the  derivation  of  proportional  gain,  proportional- 
plus-integral  and  adaptive  multiple  model  controller 
designs.  Chapter  4  discusses  the  Monte  Carlo  simulation  and 
the  analyses  to  be  performed,  which  are  in  accordance  with 
the  objectives  stated  in  Section  1.2.  Chapter  5  presents 
the  results  from  the  Monte  Carlo  analyses  along  with  a 
detailed  interpretation  of  the  findings.  One  finding,  an 
apparent  instability  problem  found  during  the  robustness 
analysis  of  the  beam  time  constant  is  developed  more  fully 
in  a  separate  chapter,  Chapter  6.  This  chapter  analyzes  the 
stability  and  robustness  characteristics  of  the 
deterministic  and  stochastic  PI  controller  designs,  the 
Kalman  and  Snyder-Pishman  filters,  and  the  Meer  filter 
structure.  The  conclusions  and  recommendations  for  future 
research  are  provided  in  Chapter  7 . 


The  first  challenge  in  designing  a  controller  for  the 
pointing  and  tracking  of  a  neutral  particle  beam  is 
developing  an  estimator  which  can  predict  the  location  of 
the  beam.  Previous  research  resulted  in  defining  the  beam 
model  and  creating  the  Meer  filter,  a  Multiple  Model 
Adaptive  Estimator  (MMAE)  composed  of  a  bank  of  Snyder- 
Fishman  filters.  This  chapter  will  explain  the  particle 
beam  model,  the  Poisson  space-time  point  process,  and  the 
Meer  filter.  The  chapter  concludes  with  a  proposal  to 
simplify  the  Meer  filter  under  a  given  set  of  conditions. 

2.1  Particle  Beam  Model 

One  recommended  method  for  locating  the  center  of  a 
particle  beam  entails  illuminating  the  beam  by  one  or  more 
lasers  [23] .  As  a  result,  the  particle  electrons  absorb  the 
photons  from  the  laser  and  jump  to  a  higher,  unstable  energy 
state,  and  then  spontaneously  decay,  returning  to  the  ground 
energy  state  and  dissipating  the  energy  as  photons.  The 
photons  radiate  approximately  in  an  isotropic  manner  and  the 
position  of  the  beam  can  be  inferred  from  the  location  at 
which  the  photons  strike  a  photo-detector  array.  The 
photoelectric  events  occur  as  a  discrete  process  at  random 
intervals.  Because  the  optical  sensors  suffers  from  dark 


currents  and  background  light,  the  sensors  will  corrupt  the 
signal-generated  data  with  erroneous,  noise-induced  events. 

Meer  noted  that  the  signal-induced  event  could  be 
modeled  as  a  Poisson  space-time  point  process  on 
[to,-)  x  RM ,  and  that  this  is  a  model  upon  which  a  Snyder- 
Fishman  filter  could  be  based. .  Each  event  has  associated 
with  it  the  time  of  occurrence  t  €  [to , ~) ,  and  a  spatial 
location  r  S  RM  .  A  physical  detector  array  will  result  in 
quantizing  of  RM  into  a  finite  number  of  sections.  For 
simplicity,  the  quantization  is  ignored  without  loss  of 
generality,  and  the  array  is  assumed  to  be  continuous;  that 
is,  it  is  assumed  that  the  value  of  jr  is  known  as  an  exact 
value  from  a  continuous  range  of  possible  values. 

The  rate  of  occurrence  is  defined  as  the  signal  rate 
parameter,  Xs  (t ,£,2c(t) )  ,  and  is  assumed  to  have  a  spatial 
Gaussian  function  given  as  [18] : 

As  (t,£,£(t)  )  *  A  ( t )  exp  { -  [£-H  ( t )  jt  ( t )  ] T  R-  1  (t)  [  r-g.  { t )  x  ( t )  ]  /  2 } 

(2-1 

where 

A ( t )  is  the  maximum  amplitude  of  the  rate  function 
£  is  a  symmetrical  positive  definite  matrix 
defining  the  spread  of  the  beam 
£  is  a  spatial  location  on  the  detector  array 
H<t)a(t)  is  the  signal-inferred  beam  centroid  in  RH 
&(t)  is  a  stochastic  process  defining  the  state 
dynamics  of  the  beam  centroid 
g(t)  is  an  m  x  n  projection  matrix  from  the  state 
space  into  the  space  of  measured  photoelectron 
events 

The  beam's  spatial  Gaussian  function  is  a  result  of  the 
beam's  diffusion. 


The  state  process  is  modeled  as  the  Gauss-Markov  output 
of  the  linear  differential  equation: 


d£{t)  -  £{t)s(t)dt  +  £(t)dfi.(t) 

S(to  )  ■  £o 


(2-2) 


where  £.(t)  is  a  Wiener  process,  whose  hypothetical 
derivative  is  white  Gaussian  noise  of  unit  strength,  and  &> 
is  a  Gaussian  random  vector  with  mean  and  covariance  £o  . 

Similar  to  the  signal-induced  events,  the  noise-induced 
events  are  modeled  as  a  Poisson  space-time  point  process  on 
[to ,•)  x  RM  with  a  noise  rate  parameter,  Xn (t,£) .  The  noise 
process  is  assumed  to  be  statistically  independent  of  the 
signal  process  [18] .  The  noise  events  are  assumed  to  be 
uniformly  distributed  over  the  detector's  field  of  view. 

The  relationship  between  the  signal-  and  noise-induced 
events  is  the  signal-to-noise  ratio,  which  is  defined  as  the 
ratio  of  the  average  signal  rate  to  the  average  noise  rate. 
Por  a  comparison  of  the  two  rate  functions,  see  Figure  2-1. 
The  signal-to-noise  ratio  is  the  ratio  of  the  areas  under 
the  two  curves  in  that  figure. 

Because  both  the  signal  and  noise  processes  are  Poisson 
point  processes  and  independent  of  one  another,  their  sum 
remains  a  Poisson  point  process.  This  can  be  shown  by 
taking  the  product  of  the  characteristic  functions  of  the 
signal  and  noise  processes.  The  resulting  characteristic 
function  will  be  a  Poisson  point  process  with  the  total  rate 
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Figure  2-1.  Signal  and  Noise  Rate  Parameter  Functions 


parameter  equal  to  the  sum  of  the  signal  and  noise  rate 
parameters: 


A(t,£.(t)  ,&(t) )  "  As  (t,£(t)  ,x(t) )  +  Am  (t,r(t)  ,jc(t) ) 


(2-3) 


>V* 
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Snyder  and  Fishman  developed  an  estimator  which 
specifically  handled  the  case  of  measurements  that  appear  as 
a  Poisson  point  process  [25] .  Their  filter  is  similar  in 
structure  to  the  Kalman  filter,  but  the  Snyder-Fishman 
filter  differs  in  two  significant  ways.  The  development  of 
the  Snyder-Fishman  filter  was  based  on  the  assumption  that 
all  measurements  were  signal-induced.  In  other  words,  the 
filter  is  limited  to  noise-free  environments.  The  other 
difference  is  that  the  sample  period  is  not  fixed  and  the 
time  between  signal  events  is  a  Poisson  distributed  process. 

With  the  absence  of  noise-induced  events,  the  Snyder- 
Fishman  filter  estimates  the  beam's  centroid  as  H(t)£(t)  , 
where  £{t)  is  the  expected  value  of  the  beam  states,  x(t) , 
conditioned  on  all  the  previous  measurements,.  The  filter 
is  described  by  the  following  differential  equations  [25] : 


.  F 


d*(t)  »  F(t)*(t)dt  +  £(t) [r-H(t)£(t)] -N(dt  x  dr) 

R" 


+  B (t)u ( t)dt  (2-4) 


d£(t)  -  £(t)£(t)dt  +  £(t)£Mt)dt  +  £(t)£Mt)dt 

£<t)H(t)P(t)N(dt  x  dr)  (2-5) 


* 

J 

Jrm 


&U)  -  E(t)iiMt)  [H<t)£(t)in  (t)  +  R(t)]-i 


(2-6) 
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where 

£(to)  ■  jco  and  g(to)  ■  go  are  the  initial  conditions 
g(t)  is  the  filter-computed  error  covariance 
£(t)£Mt)  -  g(t)fl(t)£»  (t) ;  £(t)  =  i 

£(t)  is  the  filter  gain 

The  notation  found  in  Equations  (2-4)  and  (2-5)  involves 
counting  integrals,  where 


f(t,£)N(dt  x  dr)  =  - 
JRM  It 


Nt  «  0 


(2-7) 


I  f  (ti  ,£i  )  5  (t ,  ti  )  Ht  *  1 

l-t 


and  5(t,ti)  is  a  Kronecker  delta.  Simply  stated,  if  no 
events  are  detected,  then  the  integral  equals  zero.  If  an 
event  is  detected,  the  integral  causes  a  jump  discontinuity, 
f(t,r),  to  be  added  to  the  solution  of  the  differential 
equation  for  time  ti  [25] . 

The  propagation  and  update  equations  can  be  derived 
from  Equations  (2-4)  and  (2-5).  The  propagating  equations 
for  use  between  (signal-induced)  events  are: 


d£ft)  *  F(t)£(t)dt  +  B(t)u(t)dt 


dg(t)  *  F(t)P(t)dt  +  g(t)FT  (t)dt  +  G(t)GMt)dt 


(2-8) 


(2-9) 


Equations  (2-8)  and  (2-9)  are  actually  implemented  in  the 
discrete  controller  designs  (to  be  developed  in  Chapter  3) 
as  stochastic  difference  equations: 


8<tt«l~)  -  £(tl  ♦  l  ,  tl  )«(tl  ♦  )  +  Bd  ( ti  )  u  ( ti  ) 


(2-10) 


£(ti*f)  *  £(ti«i  ,ti  )£(ti*)fiT  (ti«i  ,ti  ) 

+  Gd  <ti  )Qt  (ti  )GdMt»  )  (2-11) 

where  Gd (ti )  ■  1  and  the  calculated  as 

'ti*i 

On  (ti  )  *  •(ti*t  ,t)G(t)Gt  (T)eT  (ti  ♦  I  ,T)dT  (2-12) 

.  tl 

When  an  event  has  been  detected,  a  measurement  update  can 
take  place  as  defined  by  the  following  equations: 

S(t»M  »  k(ti-)  +  K(ti  )  [£t  -H(ti  )£(ti- )  ]  (2-13) 

£(tiM  -  P(ti-)  -  E(ti  )£(ti  )P(ti“)  (2-14) 

where  the  Snyder-Fishman  filter  gain  is  defined  by 

K(ti)  *  £(ti- )|T  (tt  )  [£(ti  )£(ti- )aT  (ti  )  +£(ti)]-»  (2-15) 

Equations  (2-8)  through  (2-15)  appear  strikingly 
similar  to  the  Kalman  filter  equations,  except  the 
measurement  update  times  are  not  known  a  priori  [25] . 
Instead,  the  measurement  updates  occur  whenever  an  event  is 
detected,  and  the  time  interval  is  defined  as  a  random 
process . 
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2.3  Multiple  Model  Adaptive 


Estimator 

The  problem  with  using  the  Snyder-Fishman  filter  is 
that  all  the  events  are  assumed  to  be  signal-induced.  To 
incorporate  the  noise-induced  events,  Meer  developed  a  bank 
of  Snyder-Fishman  filters  into  a  MMAE ,  where  each  filter  in 
the  bank  is  associated  with  an  individual  hypothesis.  The 
hypothesis  defines  which  events  are  signal-induced  and  which 
are  noise-induced.  The  result  is  a  hypothesis  tree  with 
each  branch  depicting  a  hypothesized  noise/signal  sequence, 
and  each  branch  has  an  associated  probability  that  its 
assumed  sequence  of  events  in  fact  occurred  (see  Figure 
2-2) .  In  other  words,  associated  with  each  of  these 
hypothesis  sequences  could  be  a  specific  Snyder-Fishman 
filter  that  performs  measurement  updates  when  it  receives  a 
hypothesized,  signal-induced  event  and  ignores  the 
hypothesized  noise-induced  events  (see  Figure  2-3) . 

The  state  estimate  out  of  each  filter  is  expressed  as 

&  (t)  -  E(s(t) lhj"»  ,Z»* I  (2-16) 

where  &  (t)  is  the  expected  value  of  the  beam  states,  &  (t) , 
conditioned  on  the  j-th  assumed  hypothesis  time  histories, 
hj"*,  and  the  observations  history,  Z"  ‘  ,  of  events  (ti  ,ri  ) , 
where  i  «*  1,  2,...,N».  Therefore,  the  overall  state 
estimate  of  the  MMAE  is  the  probabilistically  weighted 
average  of  the  individual  filter  state  estimates  as 
expressed  by 


Figure  2-3.  Multiple  Model  Adaptive  Estimator 


with  an  MMAE  computed  error  covariance  defined  as 


2"*-l 

P(t)  -  I  Pr[hj"»  IZ"»] -f£j  (t)  +  [«j  (t)-g(t)]  •  [Sj  ( t )  -S  ( t )  ] T  I 

j-1 

(2-1 

The  weighting  probabilities,  Pr [hj " *  I Z" 1 ] ,  are  the 
conditional  probabilities  that  a  hypothesis  sequence  is 
correct,  conditioned  of  the  measurement  history  that  has 
been  observed. 

The  weighting  probability  starts  at  to ,  with 
Pr[hj°IZ°]  »  1.  The  subsequent  weighting  probabilities 
appear  as: 


Pr  [hj  » 1  I Z" 1  ] 


is  ( tk  ,  £k  ,  8  ( tk  )  )  or  i*  ( tk  ,  £k  ) 
K  ( tk  ,  x>  ,  8  ( tk  ) ) 


•  Pr  [hj  "  *  " 


I ZM  1  ■ 1  ] 
(2-19) 


where  is  ( tk  ,£k 'S(tk  ) )  ,  in  (tk  ,£k  )  and  X  ( tk  ,£k  ,8( tk  ) )  = 
is  (tk  ,£k  ,8(tk  ) )  +  in  ( tk  , £k  )  are  estimates  of  the  signal, 
noise  and  total  rate  parameters,  and  Pr[hjNtIZNt]  is  the 
probability  of  hjMt  conditioned  on  the  most  recent  event 
that  occurred.  The  upward  branches  of  the  hypothesis  tree 
are  based  on  the  assumption  that  the  most  recent  event  was 
signal-induced,  and  the  downward  branches  of  the  hypothesis 


tree  are  based  on  the  assumptions  that  the  event  was  noise- 
induced  (see  Figure  2-2) .  The  estimate  of  the  signal  rate 
parameter,  Xs { tk ,£* ,£(t* ) ) ,  appears  in  the  numerator  for  a 
hypothetical  signal-induced  event,  and  is  found  by 
evaluating  Equation  (2-1)  for  x(t)  =  &(t) .  The  estimate  of 
the  noise  rate  parameter,  Xn{tk,£k),  is  used  for  a 
hypothetical  noise-induced  event,  and  Xn  =  Xk .  Because  of 
the  recursive  nature  of  Equation  (2-19) ,  only  knowledge  of 
the  most  recent  measurement  is  needed,  rather  than  a  growing 
memory  of  the  entire  history  of  measurements,  for  its 
evaluation  [18] . 

2.4  Pruning  the  Hypothesis  Tree 

At  this  point,  all  the  equations  of  the  Meer  filter 
have  been  derived  in  order  to  construct  a  full-scale 
particle  beam  estimator.  The  only  problem  is  that  we  have 
developed  a  growing  and  soon  to  be  unmanageable  hypothesis 
tree,  because  the  tree  must  grow  two  branches  for  each  event 
to  represent  the  possibility  of  either  a  signal  or  noise 
event.  Therefore,  if  we  have  a  total  of  Nt  events,  and  Nt 
update  cycles,  there  will  be  2N *  possible  signal-versus- 
noise  hypothesis  sequences  (frequently  referred  to  as  tree 
branches) .  The  solution  is  to  prune  the  hypothesis  tree  and 
keep  the  number  of  branches  at  a  manageable  size.  The  two 
proposed  pruning  algorithms  are  the  "best  half"  and  "merge” 


methods . 


2.4.1  II ie  Best  Half  Method  [18]  The  "best  half" 
method  was  the  method  Meer  used  in  his  dissertation  to 
achieve  an  implementable  filter.  Initially,  the  hypothesis 
tree  is  allowed  to  grow  until  it  reaches  a  prescribed  memory 
depth  of  D.  This  produces  a  tree  with  2°  hypothesis 
branches,  with  half  of  the  tree  originating  from  an  assumed 
signal-induced  event,  and  the  other  half  emanating  from  a 
noise-induced  event.  With  the  (D+l)-th  event,  the 
conditional  probability  weighting  factors  associated  with 
each  branch  are  summed  for  each  of  these  halves,  and  the 
more  probable  half  is  accepted,  and  the  branches  of  the  less 
probable  half  are  eliminated  (see  Figure  2-4) .  The 
probability  weight  is  normalized  for  the  retained  half  of 
the  tree  so  the  sum  of  the  retained  probabilities  will  equal 
one. 

The  estimates  are  propagated  forward  until  the  next 
event,  and  the  measurement  update  and  pruning  process  occurs 
again.  By  doing  so,  the  tree  never  grows  past  the  memory 
depth  of  D. 

2.4.2  The  Merge  Method  An  alternative  algorithm  to 
the  "best  half"  method  was  proposed  by  Weiss  and  associates 
[26]  and  is  designed  to  preserve  some  of  the  information 
that  would  be  lost  if  the  less  probable  half  were 
eliminated.  It  is  referred  to  as  the  "merge”  method  because 
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paired  up.  That  is,  each  pair  of  hypothesis  sequences 
differ  by  the  assumption  made  about  the  oldest  event  in  the 
current  sequence,  and  when  the  sequences  are  merged,  only 
the  oldest  event  is  dropped.  For  example,  if  the  hypothesis 
sequence  is  defined  as  hj N x  ,  where  j  =  0 , 1 , 2 ,  .  .  .  ,  2M  * -1 ,  and 
Nt  is  the  total  number  of  observed  events,  one  entire 
sequence  could  be  written  as 


It  It  It 

hj  -  (hj  (1)  ,  hj  (2) 


It 

hj  (Nt  )  I 


=  (1,0 . II 


(2-20) 


All  signal-induced  events  are  represented  by  a  "1"  while 
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noise-induced  events  are  labeled  with  a  "0”.  Continuing 
with  the  example,  it  will  be  assumed  that  the  memory  depth 
is  limited  to  two  (D=2) ,  and  the  total  number  of  observed 
events  is  two  (Nt=2).  Then  the  hypothesis  sequences  would 
be: 


< 

ha2 

=  {ha2  (1) ,  ha2  (2) ) 

=  {1,1} 

V 

ha  2 

=  {ha2  (1) ,  ha2  (2) } 

=  {1,0} 

1  * 

hi2 

=  {hi*  (1) ,  hi2  (2) } 

r) 

O 

II 

ho2 

*  {ho2  (1) ,  ho2  (2) } 

*  {0,0} 

(2-21) 

As  the  next  event  is  observed,  Nt=3,  and  there  would  be 
eight  hypothesis  sequences.  If  they  are  paired  up  by  the 
D-l  most  recent  events,  the  sequenced  pairs  would  be 


Pair 

1: 

hr3 

* 

{  [ha2 

(1),  ha2 

(2)]  , 

h7  3 

(3)  } 

= 

{1,1,11 

1 

ha  3 
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{  [hi2 

(1)  ,  hi2 

(2)]  , 

ha3 

(3)  } 

= 

{0,1,1} 

Pair 
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he3 
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{  [ha2 
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ha3 

= 
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(2)]  , 

ha3 
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(3)  } 
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{0,1,0} 
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h« 3 

(3)  } 
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{1,0,0} 

fed 

ho3 

s> 

{  [ho2 

( 1 )  ,  ho 2 

(2)1  , 

ho3 

(3)  I 

s 

{0,0,0} 

(2-22) 


To  merge  the  paired  hypothesis  sequences,  the  oldest,  i.e, 
(Nt-D)-th,  event  is  dropped,  and  the  new  weighted 
probability  is  the  sum  of  the  probabilities  of  the 


24 


individual  sequences  that  were  merged.  The  state  estimates 
and  error  covariances  associated  with  each  merged  pair  can 
be  found  using  the  following  equations: 


’ (t)  ■  {Pr[hj"‘  |ZMt  ]  •&  <t)  +  Prthk**  IZNt  ]  -gk  (t)  } /A  (2-23) 

Pj*(t)  -  {Pr[hj»MZ"»]  •  (£j  (t)  +  [&  (t)-gj ' (t)  ]  [&i  (t)-Si  '  (t)]M 
+  Prthk"*  !Z«*]  •  (Pk  ( t)  +  [gk  (t)-gj  *  (t)]  [gk  (t)-gj  •  (t)]M  I/A 

(2-24) 

A  *  Pr[hj  |Z“*  ]  =  Pr[hj»*  !ZNt  ]  +  Pr  [hkM »  |  Z"  *  ]  (2-25) 


where 


hj D  and  hkD  denote  the  two  different  sequences 

within  each  pair  (see  Equation  (2-19) :  k  *  j  +  2° 
and  j  =  0,1,.., (2° -1) 

D  is  the  memory  depth  of  the  hypothesis  tree  after  the 
pruning  process  has  occurred 

'  and  £4  '  are  the  "merged"  state  estimates  and 
error  covariances  as  the  number  of  elemental  filters 
are  reduced  from  20+l  to  2® ;  j'  =  0 , 1 , . . . , (2® -1) 

Equation  (2-25)  must  be  normalized  so  that  the  sum  of  the 

probabilities  of  each  new  estimate  equals  one.  Figure  2-5 

demonstrates  the  "merge"  method.  This  method  conceptually 

has  an  advantage  over  the  "best  half"  method  because  it 

accounts  for  all  possible  time  histories  rather  than 

deleting  half  of  the  branches  from  a  decision  tree. 

Unfortunately,  it  is  computationally  more  burdensome. 


Figure  2-5. 
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2.5  SvnpUfiy.lhg  the  Mgey  Filter 

Zicker  ran  a  performance  analysis  on  the  Meer  filter 
and  found  the  Meer  filter  virtually  insensitive  to  the  depth 
parameter  D.  Varying  D  from  D*1  to  D*8  produced  less 
than  one  percent  change  in  rms  errors.  This  is  partially  a 
function  of  using  simple  scalar  first  order  Gauss-Markov 
model  to  define  the  beam  centroid  dynamics  [17]  .  Assuming 
the  first  order  beam  model  proves  to  be  adequate,  this 
insensitivity  could  be  used  to  simplify  the  MMAE  structure 
of  the  Meer  filter  by  limiting  the  depth  to  D*l,  and 
defining  the  filter  gain,  &(£.*  )  ,  as  a  function  of  the 


*  v 


!  U 


The  probability  it  was  in  fact  a  signal-induced  event  is 

Pi  (ti  )  -  Xs  (ti  ,£i  ,£<ti  )  )  /  X  (ti  ,£i  ,&(  ti  ) )  (2-27) 

using  the  notation  of  Equation  (2-19) .  If  the  elemental 
Snyder-Fishman  filter  assumed  that  the  event  was  noise- 
induced,  then  &a(tiM  =  £a(tf).  The  probability  it  was  a 
noise-induced  event  is 

pa  (ti  )  =  Xn  (ti  ,ri  )  /  X  (ti  ,£t  ,£(ti  ) )  (2-28) 

Therefore,  the  adaptive  beam  state  estimate  is 

S*(tiM  *  pi  (ti  )*i  (ti  *  )  +  pt  (ti  )£s  (ti*  )  (2-29) 

-  Sa  (ti-  )  +  pi  (ti  )£(ti  )  [£  -  Hi  (ti  )Xi  (tA-  )  ] 

where  &A(ti-)  *  3jti(tf)  «  Sa(ti~).  Thus,  the  algorithm 
involves  a  single  Snyder-Fishman-like  filter,  but  with  a 
modified  gain  of  [pi  (ti)K(ti)],  and  thus  the  gain  is  a 
function  of  the  residual  as  seen  in  Figure  2-6. 

Higher  order  Gauss-Markov  beam  models  should  have  a 
greater  sensitivity  to  D.  Therefore,  this  simplification 
may  be  limited  to  this  simple  case. 


2.7  Summary 
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The  Meer  filter  is  a  Multiple  Model  Adaptive  Estimator 
(MMAE)  that  is  based  on  a  space-time  point  process  model  for 
measurement  events,  and  it  can  be  used  to  estimate  the 
position  of  a  neutral  particle  beam.  The  Meer  filter 
locates  the  beam  by  detecting  the  photo-electric  event 
occurring  from  the  spontaneously  decaying  excited  particle 
within  the  beam.  The  time  between  events  is  a  Poisson 
process,  and  this  is  properly  modeled  by  the  Snyder-Fishman 
filter.  Meer  introduced  the  MMAE  structure  to  handle  the 
noise-induced  events  which  the  Snyder-Fishman  filter 
neglected.  Because  there  is  no  deterministic  mechanism  for 
declaring  whether  an  event  is  signal-  or  noise-induced,  a 
hypothesis  tree  was  developed.  The  depth  of  the  tree  was 
limited  by  incorporating  either  the  "best  half"  or  "merge" 
method.  The  results  from  Zicker's  research  indicated  that 
the  Meer  filter  could  be  reduced  by  limiting  the  depth  to 
one,  redefining  the  filter  gain  as  a  function  of  the 
residual,  and  eliminating  the  multiple  model  structure. 

With  the  beam  state  estimator  in  hand,  the  next  step  is 
to  develop  the  controller (s)  to  regulate  the  beam  and  track 
the  target.  The  next  chapter  will  address  this. 


The  next  challenge  is  to  design  a  controller  which  can 
place  the  beam  on  the  target.  Zicker  developed  the  first 
controller  for  this  application  [27]  .  It  was  a  proportional 
gain  controller  and  was  designed  as  an  intuitive  tool  to 
gain  insight  into  the  tracker  problem.  But,  because  the 
proportional  gain  controller  does  not  have  type-one 
characteristics  that  allow  rejection  of  unknown,  constant 
disturbances  that  arise  from  linearized  models.  Moose 
developed  a  proportional-plus-integral  (PI)  controller 
[20].  Jamerson  carried  on  Moose's  research,  and  developed  a 
different  PI  controller  based  on  a  more  realistic  target 
model  [8] .  Although  the  PI  controller  could  reject  constant 
unmodeled  disturbances,  it  did  so  at  the  cost  of  additional 
sluggishness  and  had  difficulty  tracking  a  highly 
maneuverable  target.  The  next  step  will  be  to  develop  a 
Multiple  Model  Adaptive  Controller  (MMAC)  to  handle 
parameter  uncertainties  caused  by  a  highly  manueverable 
target  that  can  display  rapid  variations  from  a  benign 
trajectory  to  evasive  maneuvering. 

All  the  controllers  are  based  on  the  same  fundamental 
assumptions:  a  linear  system  model,  a  deterministic  optimal 

control  law,  full  state  feedback,  and  a  quadratic  cost 
criterion.  The  regulator's  gain  is  based  on  a  backward 
Riccati  difference  equation,  assuming  all  the  particle  beam 


states  are  perfectly  known.  Then,  by  assumed  certainty 
equivalence,  the  full  state  feedback  is  replaced  by  a 
filtered  or  "observed"  state  feedback  provided  by  the  Meer 
filter.  The  tracker  portion  of  the  controller  is  designed 
around  the  regulator,  and  directs  the  beam  from  the 
regulator's  zero  setpoint  onto  the  target,  using  target 
state  estimates  generated  by  a  Kalman  filter. 

Because  the  target  model  is  used  to  developed  the 
trackers,  it  will  be  presented  first,  followed  by  the 
proportional  gain,  proportional-plus-integral  and  the 
multiple  model  adaptive  controller  designs.  The  controllers 
are  developed  in  the  general  vector  form  and  then  reduced  to 
the  specific  one-dimensional  tracker  problem  model  that  will 
be  used  during  the  analyses. 

3.1 

Zicker  [27]  and  Moose  [20]  designed  their  controllers 
for  a  simple  target  model  composed  of  a  first  order  Gauss- 
Markov  position  process.  Research  centered  around 
feasibility  studies  and  the  emphasis  was  on  gaining  the 
knowledge  and  insight  required  to  develop  a  realistic 
controller.  Starting  with  Jamerson  [8] ,  the  emphasis 
shifted  to  a  more  realistic  target  model  built  around  a 
first  order  Gauss-Markov  acceleration  process.  The  linear, 
time-invariant  state-space  representation  of  this  model  is 


*rr  (t) 
(t) 
*t  a  ( t ) 


0  10 
0  0  1 

0  0  -1/tt 


Xtp  (t) 
XT  *  (t) 
XT  A  (  t  ) 


WT  (t) 


(3-1) 


where 


xt  r  (t)  is  the  target  position  state 
xt v (t)  is  the  target  velocity  state 
Xt  a ( t )  is  the  target  acceleration  state 
tt  is  the  acceleration  correlation  time  constant 
wt (t)  is  a  zero  mean  white  Gaussian  noise  of  strength 
Qt  associated  with  the  target 
£  is  the  matrix  that  describes  the  system's  dynamics 
£  is  the  matrix  that  maps  the  white  noise  effects  into 
the  state  vector. 


Equation  (3-1)  has  three  poles,  of  which  two  are  at  the 


origin  of  the  s-plane.  This  equation  of  target  motion  is 


inherently  astable,  and  the  target  is  uncontrollable  by  the 


control  input,  u(t);  therefore,  the  steady  state  Riccati 


equation  used  to  generate  controller  gains  may  not  have  a 


solution  [10,14].  All  modes  must  be  stabilizable  to 


guarantee  a  solution  of  the  steady-state  Riccati  equation. 


To  insure  a  solution,  the  following  sections  will  move  the 


two  poles  at  the  origin  of  the  s-plane  to  the  left  by  some 


small  epsilon.  This  will  guarantee  a  solution  for  the 


steady  state  Riccati  equation,  but  not  affect  the  filter's 


development . 


In  order  to  use  the  target  equation  of  motion  defined 


by  Equation  (3-1) ,  the  target  position  will  be  measured  as  a 


discrete  function  defined  by 


z(ti  )  ■  fr  (ti  )gr  (ti  )  +  v(ti  ) 


■  -  'MVA'.sW 


SI 


XT  P  (  ti  ) 

z(til  »[100]  xtv(ti)  +  v(ti) 

XT  A  ( tl  ) 


(3-2) 


where  [  xtp  (ti  )  xtt  (ti  )  XT*(ti)  ]T  is  the  target  state  and 
v(ti)  is  a  discrete-time  zero-mean  measurement  corruption 
noise  with  a  covariance  R,  and  is  assumed  to  be  independent 
of  the  dynamics  driving  noise  in  Equation  (3-1) .  How  the 
target  is  actually  detected  and  located  is  beyond  the  scope 
of  this  thesis.  Between  measurements,  the  state  estimate 
and  covariance  matrix  are  propagated  forward  in  time  by 


St  ( ti  ♦  1  “  )  *  0T  ( ti « 1  ,  ti  )  £t  ( ti  *  ) 


£t  ( ti  ♦  i  -  )  ■  4>t  ( ti  ♦  t  ,  ti  )  £t  ( ti  ♦  i  )  <&t  T  ( ti  ♦  i  ,  ti  ) 


(3-3) 
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Or  (ti  ♦  i  ,  t)Gt  (t)2t  (t)£tt  (t)Ott  (ti  ♦  i  ,  t)dt  (3-4) 


where  £r(tt*i,ti)  is  the  state  transition  matrix  associated 
with  £  in  Equation  (3-1) .  During  the  measurement  update, 
the  state  estimate  and  covariance  matrix  are  updated  by 


St  (ti*  )  =  St  (ti-  )  +  £(ti  )  [Z(ti  )  -  Ht  (ti  )St  (tf)]  (3-5) 


EtUiM  -Pr(ti-)  -  £(ti  )&t  (ti  )Pt  (ti- )  (3-6) 


f£(ti  )  »  Pt  (tf  )|!tt  (ti  )  [Hr  (ti  )PT  (ti- )HrT  (ti  )  +RT(ti)]-» 
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where  the  initial  conditions  were 
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and  the  subscript  "T"  denotes  "target." 


As  with  any  Kalman  filter,  if  the  measurement  noise 


variance  R(ti)  is  increased,  then  the  filter  must  rely  more 


on  the  dynamics  and  propagation  equations.  In  a  similar 


fashion,  if  the  driving  noise  strength  Qt  is  increased,  the 


filter  must  rely  more  heavily  on  the  incoming  measurement. 


Gain  Controller 


If  we  had  a  linear  quadratic  Gaussian  (LQG)  stochastic 


controller  problem,  we  could  use  the  certainty  equivalence 


property.  But,  the  adapative  mechanism  in  the  Meer  filter 


violates  linearity;  therefore,  we  must  use  the  assumed 


certainty  equivalence  methodology  [14] .  It  allows  us  to 


synthesze  a  sub-optimal,  non-linear  stochastic  controller  by 


generating  the  associated  optimal  deterministic  full-state 


feedback  controller,  and  then  replacing  the  actual  states 


with  the  conditional  mean  estimates  from  the  stochastic 


filters . 


3.2.1 


Gain 


The  function  of  the 


regulator  is  to  drive  the  states  of  the  beam  to  a  zero 


setpoint  defined  as  the  center  of  the  detector  array.  The 


linear  discrete-time  state  equation  that  defines  the  beam  is 


&  ( ti  ♦  i  )  *  ( ti  + 1  ,  ti  )&  ( ti  )  +  &d  ( ti  )up  ( ti  )  +  Gd  w«t  (ti  )  (3-8) 


where 


■J 

*  « * 


s  i 


&■  (ti  )  is  the  position  of  the  bean 

u(ti )  is  the  control  applied  over  the  next  period  from 
ti  to  ti  ♦  1 

Wd (ti )  is  a  zero  mean  white  Gaussian  discrete-time 
stochastic  process  of  strength  Qd ,  where 


ti  ♦  i 


$b  ( ti  ♦  1  ,  T )  Gb  {  T  )  Qfi  ( T )  Gb  T  ( T )  Ob  T  (  ti  ♦  1  ,  T  )  dT 


(Note  that  the  Snyder  and  Fishman  filter  defines  dynamics 
driving  noise  strength  as  G(t)GT(t);  therefore,  Qp(t)=I.) 
The  discrete,  optimal  deterministic  control  law,  assuming 
perfect  knowledge  of  jcb  (ti  )  ,  is 


H*  (ti  )  «  -Gc*  (ti  )jcb  ( ti  ) 


(3-9) 


where  Gc* (ti )  is  the  optimal  controller  gain.  Because  we  do 
not  have  perfect  knowledge  of  the  beam  position  state, 

2Cb  (ti  )  ,  it  is  replaced  by  the  conditional  mean  of  the  beam 
position  state,  ka (ti  ) ,  which  is  provided  by  the  Meer 
filter.  This  is  done  in  accordance  with  the  assumed 
certainty  equivalence  synthesis  methodology.  The  result  is 
a  discrete,  optimal  stochastic  control  law: 


H*  (t»  )  *  -£c*  (ti  )£b  (ti  ) 


(3-10) 


The  controller's  optimality  is  defined  by  minimizing  a 
quadratic  cost  function  defined  as 


■'ts- 


J  -  Ef  Z  H[&T  <ti  )&<ti  ) 2£( ti  )  +  uMti  )U(ti  )u(ti  )] 

4*0 


+  J4xT  ( tM  ♦  i  )  JJf  ( t«  ♦  i  )  £ ( t*  ♦  i  )  ( 3-11 ) 


where 


&(tN«i)  is  the  final  position  to  be  achieved 
$(ti )  is  the  position  at  time  ti 
&  and  Xr  are  n-by-n  constant  weighting  factors 

defined  as  a  positive  semi-definite  matrices 
U  is  an  r-by-r  weighting  factor  defined  as  a  positive 
definite  matrix 

The  controller  gain,  Gc(ti),  in  Equation  (3-10)  is  generated 
by  solving  the  backwards  Riccati  recursion 


GcMti)  »  [U(ti  )  +  &dT  (ti  )RC  (ti*i  )Bd  (ti  )]-* 

•  [gd T  ( ti ) ».  c  ( ti  4 1 )  ( ti  ♦  i  ,  ti )  ]  (3-12) 

Kc(ti)  »£(ti)  +  eMti  4  I  ,  tt  )Kc  (ti  4  1  ) 

•  [^T  (t!4i  ,ti  )  -  Id  ( ti  )  ]  .Gc*  (ti  )  (3-13) 

solved  backwards  from  the  terminal  condition 


Kc  (t*4i  )  =  Xf  (3-14) 

The  resulting  equations  are  sufficient  to  define  the  beam 
regulator  as  depicted  in  Figure  3-1. 

For  simplicity,  the  problem  will  be  reduced  to  one 
physical  dimension.  That  is,  the  beam  and  target  models 
will  be  confined  to  one-space,  R1 ,  and  Bd ,  Ga  and  X(ti ) 
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Figure  3-1.  Proportional  Gain  Regulator 

will  be  scalars  equal  to  one.  The  reason  the  weighting 
matrix  X(ti  )  can  be  reduced  to  unity  is  that;  the  X/U  ratio 
is  the  important  factor  in  determining  the  steady  state 
gains,  and,  for  the  scalar  case,  either  the  numerator  or 
denominator  can  be  set  to  one,  and  the  other  weighting 
function  adjusted  to  maintain  the  proper  ratio.  This 
simplification  is  limited  to  the  scalar  case. 

With  these  simplifying  assumptions,  the  general 
regulator  equations  reduce  to: 


xa(ti«i)  *  0*x»  (ti  )  +  u(tj) 


(3-15) 


v *  *  •  *  v'n. V'” 


*  exp [- ( ti  ♦  i -ti  ) /tb  ]  =  exp[-At/tB] 


(3-16) 


u(ti  )  =  -Gc  *  ( ti  )  XB  ( ti  ) 


(3-17) 


Si 


Gc*  (ti  )  »  Kc  (ti*i  )<I>b  /  [U  +  Kc  (ti  +  t  )] 


(3-18) 


Because  this  is  a  pure  tracking  problem  without  a  well 


defined  terminating  point,  the  solution  in  found  by 
determining  the  steady  state  solution  to  the  backward 
Riccati  difference  equation  (i.e.,  solving  for  the  positive 
root  in  a  quadratic  solution) : 


Kc  (ti  )  =  Kc  (ti  ♦  i  )  ■  Kc 


Kc2  +  Kc  [U(l  -  $b2  )  -  1]  -  U  =  0 


(3-19) 


3.2.2  Proportional  Gain  Tracker  Now  we  wish  to 
extend  the  regulator  to  a  tracking  problem  using  the  same 
LQG  synthesis  process.  The  goal  is  to  minimize  the 


difference  between  the  controlled  beam  variables,  yc  ,  and 
the  target  or  reference  variables,  yr .  The  variables  are 
the  linear  transformations  of  the  states  of  the  beam  and 


target  states  expressed  as: 


(ti  )  *  £c  (ti  )xb  (ti  ) 


(3-20) 


2,r  ( 1 1  )  =  Cr  ( tl  )  ST  (tl  ) 


(3-21) 


$ 


and  the  difference  between  them  is  referred  to  as  the 
tracking  error,  expressed  as 


e(ti  )  =  Cc  (ti  )x»  (ti  )  -  Cr  (tl  )XT  (ti  ) 


[  Cc  (tl  )  -£r  (tl  )  ]  [  Xb  (tl  )  XT  (tl  )  ] 


Ca  ( tl  )  Xa  ( tl  ) 


(3-22) 


The  tracking  error  shall  be  regulated  to  zero  by  way  of 
minimizing  the  tracker  augmented  cost  function  defined  as 


J  »  Ef  i  Ji[x-T  {tl  )x,  (tl  )£,  <tt  }  +  ut  (tl  )u(ti  ) u ( 1 1  )] 


♦  (t*»  1  )Xfa  (tN.  I  )2Ca  (ta*i  )  (3-23) 


where 


&  (ti  )  -  £aT  (ti  )*•  ( ti  )£.  (ti  ) 

Xf  a  ( 1 1  )  *  £a  T  ( t*  ♦  l  )  JTf  •  ( t*  ♦  l  )  £a  ( ta  ♦  l  ) 


Assuming  we  will  have  perfect  knowledge  of  both  the  beam 
states  and  target  states,  we  can  obtain  the  full-state 
proportional  gain  control  law  of 


H*  (ti  )  -  -  [  gciMti  )  Gd'(ti)  ]  [  2£b  (ti  )  xt  (ti  )  ]t 

■  “  £c  1*  (tl  )£•  (tl  )  -  gc2*  (tl  )&T  (tl  ) 


(3-24) 


The  solution  to  the  augmented  gain  matrix  is  solved  by  an 
augmented  backward  Riccati  recursion: 
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(3-25) 

Ecn(ti)  -  £cT  (tl  )X(tl  )£.C  (tl  )  +  icTScn  ( t»  ♦  1  )Oc 

-  icTKcn  (tiM)£d(ti)5ci‘(ti)  (3-26) 

Ecl*(tl  )  »  -£cT  (tl  )X(tl  )£r  (t»  )  +  acTSp»*  <tl*l  )Or 

~Gcl*Mtl  )  [  £dT  (tl  )Scn  (tiM  )±r  ]  (3-27) 

where  <fec  is  the  beam  state  transition  matrix,  0»(ti»i,ti) 
is  this  application,  Or  is  the  target  state  transition 
matrix,  <fcr ( ti  ♦  i  , ti  ) ,  and  Equations  (3-26)  and  (3-27)  are 
solved  backwards  from  the  terminal  conditions 


Sen  (tut  i)  =  £cT  (tH»i  )Xf  (tN  +  i  )Cc  (ti.fi  )  (3-28) 


Ec  n  ( t»  ♦  i  )  *  -Cc  T  ( tn  ♦  i  )  Yf  ( tN  ♦  i  )  Cr  ( tN  ♦  i  )  ( 3-29 ) 
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By  carefully  partitioning  the  augmented  backward 


Riccati  recursion  equations,  one  can  see  that  the  solution 


for  GeiMti)  is  identical  to  the  feedback  gain,  G *  ( 1 1  )  , 


found  in  the  regulator  (Note:  Equations  (3-25)  and  (3-26) 


are  similar  to  Equations  (3-12)  and  (3-13)).  Thus,  Ken  and 


Gci*  are  independent  of  the  reference  variable  being 


tracked,  and  they  are  found  by  solving  the  deterministic 


regulator  design.  The  tracker  gain,  Gc 2 ,  is  determined  from 


Equations  (3-27)  and  (3-29) ,  and  is  a  function  of  the 


feedback  gain  and  the  reference  variable: 


SpiMti)  -  (U(ti)  +  Bd T  ( ti  ) Kc  1 1  ( ti  ♦  1  ) Bd  ( ti  )  I ~  1 


£d  T  ( ti  )  Kp  1  2  ( ti  ♦  1  )  <£r  ( ti  +  1  ,  ti  )  (3-30) 


The  resulting  tracker  model  is  diagrammed  in  Figure  3-2. 


Once  again,  the  thesis  problem  will  be  reduced  to  one 


dimension,  and  the  tracker  error  will  be  the  scalar 


difference  between  the  beam  and  target  positions.  As 


foretold,  the  regulator  design  will  provide  the  feedback 


gain,  Gc 1  * .  The  problem  setup  is  as  follows: 


I  »  Xu 


(3-31) 


Sfi  -  1 


(3-32) 


£r  -  [  1  0  0  ] 


(3-33) 


C.  -  [  Cc  -Cr  ]  =  [1-1  0  0  ] 


(3-34) 
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Xb  (ti  ) 

uMti)  -  -[  Gc  l  *  GctP*  Gc  *  v  *  Gcja‘  ]  Xtp(ti) 

Xt?  (ti  ) 
XT  A  (tl  ) 


(3-36) 


The  form  of  the  augmented  weighting  matrix,  X» ,  can  be 
intuitively  verified.  We  are  trying  to  drive  the  quantity 


e(ti  )  »  xb  (ti  )  -  xtp  (ti  ) 


(3-37) 


towards  zero,  so  the  quadratic  penalty  placed  on  the  error 


KXne(ti)2  =  M  xb  (ti  )  XTp(ti)  ] 


Xi  i  -Xi  i  Xb  ( ti  ) 


-  Xi  i  Xi  1  XT  p  ( ti  ) 


(3-38) 


which  is  the  upper  left  block  of  Equation  (3-35) .  The  rest 
of  the  elements  in  should  be  zero  because  we  wish  to 
disregard  the  velocity  and  acceleration  components  of  the 
target  when  we  calculate  the  error  and  the  control  input. 

By  restricting  the  tracker  design  to  a  time-invariant 
system  with  a  target  model  driven  by  a  stationary  noise  and 
by  implementing  constant  cost  matrices,  we  can  solve  for  the 
steady  state  tracker  gain  Get*.  As  in  the  regulator 


design,  B<i  is  set  to  one.  The  beam  state  transition  matrix 
is  a  scalar  function  decribed  by  Equation  (3-16) ,  where  tb 
is  the  time  constant. 

The  target  state  transition  matrix,  Or(ti  +  i,ti),  must 
be  stable  to  find  the  steady  state  tracker  gain,  since  the 
augmented  system  must  be  stabilizable .  If  the  target  state 
transition  matrix  has  any  poles  on  or  outside  the  unit 
circle  (z-plane  analysis),  then  the  augmented  system  model 
will  have  modes  that  are  astable  or  unstable,  and 
uncontrollable.  Because  the  original  target  state 
transition  matrix  would  be  astable,  we  have  to  move  the  two 
poles  at  the  origin  of  the  s-plane  left  by  some  small  e. 

The  modified  £  matrix  from  Equation  (3-1)  is 

-«  1  0 

0  -«  1 

0  0  — 1/T 

The  resulting  state  transition  matrix  is 

♦t  (ti*»  ,ti  )  «  J-i  f  t  si  -  £'  ]-»  I  -  (3-40) 

where  represents  the  inverse  Laplace  transformation. 

Previous  experience  found  that  the  value  of  c  *  .0001 
ensured  at  least  an  order  of  magnitude  difference  between 
the  poles  forced  away  from  the  origin  and  the  exponential 
time  constants  of  either  the  target's  acceleration,  tt ,  or 
the  beam's  position,  tb . 


(3-39) 


As  a  result,  the  tracker  gain  equations  can  be  reduced 


i 

to 


Ge  «  *  *  Kc  u  $r  /  (O  +  Kell  ) 


(3-41) 


£cia  m  [Xii  0  0]  +  T  Kc  i  a  -  Gc  i  *  T  |Cc  i  2  4>t  ( 3-42 ) 


3.3  Proportional -Plus -Integral  Controller 

Because  a  proportional  gain  controller  could  not  handle 
constant  unmodeled  disturbances  in  the  system.  Moose 
developed  the  type-one,  proportional-plus-integral  (PI) 
controller.  The  PI  controller  design  includes  many  of  the 
same  assumptions  used  in  the  proportional  gain  controller. 
Both  the  PI  regulator  and  tracker  are  designed  on  the  basis 
of  linear  system  models  and  quadratic  cost  functions,  and 
they  both  implement  steady  state  controller  gains.  Once 
again,  assumed  certainty  equivalence  is  used  and  the  Meer 
filter  provides  the  particle  beam  state  estimate  while  a 
Kalman  filter  provides  the  target  state  estimate. 

3.3.1  Proportional -Plus- Integral  Regulator  The 
function  of  this  regulator  is  to  guide  the  beam  state  to 
some  definable  setpoint,  yn .  This  setpoint  can  be  non-zero 
but,  for  this  application,  the  setpoint  will  be  defined  as 
zero.  To  handle  a  definable  setpoint  and  unmodeled  step 
disturbances,  the  regulator  augments  a  set  of  pseudointegral 
states  to  the  original  plant  state  equation.  Because  we  are 
using  a  discrete  system,  a  true  integration  process  is  not 


« 


s 


£ 


possible.  Instead,  a  pseudointegral  or  summation  process  is 
used  to  provide  a  type-one  control,  able  to  drive  the  steady 
state  mean  of  the  regulation  error,  [yc  ( ti  ) -yr  ( ti  ) -yd  ]  ,  to 
zero,  even  in  the  face  of  unmodelled  constant  disturbances 
affecting  the  plant.  Note  that  for  the  regulator,  the 
reference  state,  yr (ti  )  is  equal  to  the  zero  vector.  The 
pseudointegral  is  defined  as 

t  - 1 


g(ti  )  -  g(to  )  +  I  tyc  (tj  )  -  yr  (tj  )  -  ya  ] 

J  *  O 

»  g( ti  -  i  )  +  [yc  ( ti  )  -  yr  (ti  )  -  ya  ] 


(3-43) 


Figures  3-3  and  3-4  demonstrates  how  the  pseudointegral  is 
implemented  in  the  PI  controller. 


The  PI  control  law  is  found  by  augmenting  the 
pseudointegral  (Equation  (3-43))  with  the  beam  state 
(Equation  (3-8))  to  form 


&  ( tl  ♦  1  )  =  Oa  X a  ( ti  )  +  Bd  a  U  ( ti  )  +  Da  yd  +  Gd  a  Wd  ( ti  ) 


x»  ( ti ♦  i  ) 
a(ti  ♦  i ) 


Ob  0  JC b  ( t  i  )  Bd  0.  Gd 

+  U  ( ti  )  -  yd  +  Wd  (ti  ) 

Cc  I  g(ti)  0  io 


(3-44) 


The  cost  function  retains  the  same  structure  as  before,  but 
includes  a  quadratic  weighting  on  the  pseudointegral  in 
augumented  form: 


J  *  E(  I  K[x.T  (ti  ) Xa  (ti  )&,  (ti )  +  uT  (ti  )U(ti  ) u ( 1 1  )] 

1*0 

)i  [jCa  T  (  ta  ♦  1  )  JCf  a  (  tw  ♦  1  )  Xa  ( til  ♦  l  )  )  (3"45) 


The  optimal  control  law 


H(ti  )  =  — Gca*  (ti  )Xa  (ti  )  +  Eyd 


(3-46) 


minimizes  the  cost  function.  The  solution  to  the  optimal 


n 

5  fe 


£  -  t  gpiMti  )  -  GcaMti  )Kcn-Mti  )Kc»aT  (ti  )  ]flu  +  £22 


(3-48) 


where 


flu 

It  * 

[*»  (ti  *i  ,  ti  )  -  I] 

Bd 

Am 

£** 

Cc 

0 

(3-49) 


&  (ti )  is  derived  from  the  backward  Riccati  difference 
equation 

1 1  ( ti  )  £c  1  2  ( ti  ) 

£c  (ti  )  * 

gcl2T  (tl  )  £c  2  2  ( tl  ) 

**  £•  +  T  ( ti  ♦  1  ,  ti  )  £c  ( ti  ♦  1 )  ( ti  ♦  1  ,  ti  ) 

-  [£a«T£c  (ti«i  )0m  ( ti  ♦  1  ,ti  )]TGcMti  )  (3-50) 

solved  backwards  from  Kc  (t**i )  ■  Xf» . 

If  we  assume  a  time-invariant  system  model  and  use 
constant  weighting  matrices,  the  regulator  can  be  designed 
with  steady  state  controller  gains  (provided  the  gain 
transients  are  short  compared  to  the  total  time  interval  of 
interest,  which  is  the  case  for  this  application) .  Because 
the  tracker  design  will  develop  different  controller  gains, 
the  specialization  to  the  one-dimensional  problem  has  been 
omitted. 


tracker,  we  want  to  use  the  same  PI  controller  structure  as 


found  in  Figure  3-4  while  implementing  the  tracking  error 
equation  used  for  the  proportional  gain  tracker, 

§.(tl  )  *  Cb  (tl  )Xa  (ti  > 

*  [  Cc(ti)  -Cr(ti)  ]  [  xiMti)  £TT(ti)]T  (3-51) 

where  C«  is  an  augmented  matrix,  and  s>  is  an  augmented 
vector  that  include  the  beam  state  and  the  target  state. 
Unless  we  wish  to  offset  the  beam  from  the  target  centroid 
(as  to  track  some  other  appropriate  point  on  the  target  or 
to  lead  the  target),  the  set  point,  yd ,  will  be  defined  as 
the  zero  vector.  The  remaining  tracker  development  closely 
follows  the  PI  design. 

The  PI  control  law  is  determined  from  the  discretized, 
augmented  system  equation 


£•  (  tl  ♦  l  )  *  2Ci  (  1 1  )  +  Bd  ■  U  ( 1 1  )  +  £■  yd  +  Gd  a  Wd  ( 1 1  ) 


5p  (ti ♦  i  ) 

■  « 

Ob  0  0 

Si  (tl  ) 

Bd 

0 

Gd 

St  ( ti  -  1  ) 

■ 

o  Ot  2. 

St  (ti  ) 

+ 

0 

u  ( ti  )  - 

0 

£d  + 

a 

g(ti  ♦  i ) 

£c  -Cr  I 

g(ti ) 

m  m 

0 

I 

0 

(3-52) 


where 


Equation  (3-52)  is  composed  of  the  beam  state  Equation 
(3-8) ,  the  discretized  form  of  the  tracker  state  equation 
(3-1)  and  the  pseudointegral  Equation  (3-43)  in  which  yc  (ti ) 
is  expressed  in  terms  of  &(ti).  As  mentioned  before,  the  £ 
matrix  of  Equation  (3-1)  is  astable  and  the  poles  will  be 
moved  to  the  left  of  the  s-plane  origin. 

Again,  the  cost  function  retains  the  same  basic 
structure,  but  now  it  includes  both  the  tracker  and  the 
pseudointegral  cost  terms.  Because  we  are  using  the  same  PI 
control  structure  defined  by  Figure  3-4,  the  optimal  control 
law.  Equation  (3-46),  and  the  solution  to  the  optimal  gain 
functions.  Equations  (3-47)  through  (3-50),  are  still  valid 
and  will  be  used  to  design  the  PI  tracker. 

As  was  done  with  the  the  proportional  gain  tracker,  the 
PI  tracker  will  be  restricted  to  a  one-dimensional,  time- 
invariant  system  using  constant  weighting  matrices.  The 
first  restriction  simplifies  the  problem  substantially.  The 
controlled  variable,  y.c  ,  the  reference  variable,  ( t i  )  ,  and 
the  setpoint,  £d ,  are  scalars  defined  as 


yc  (ti  )  »  Xb  (ti  ) 
yr  (ti  )  -  Xtp  (ti  ) 


(3-53) 

(3-54) 


(3-55) 


The  latter  two  restrictions  along  with  a  decision  to  ignore 
short  transients  in  Riccati  solutions  allow  the  problem  to 
be  solved  using  steady-state  constant  gains. 


The  PI  tracker  problem  setup  is  as  follows  (see 


Equations  (3-20)  and  (3-21) : 


Ca  =  (  Cc  -Cr  ]  =  [  1  :  -1  0  0  ] 


(3-56) 


The  state  transition  matrices  of  the  beam  and  target  are 


<$b  *  exp[- ( ti  ♦  i  -ti  ) /tb  ]  =  exp[-At/tB] 


(3-57) 
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The  augmented  transition  matrix  is 
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To  be  consistent  with  Jamerson’s  derivation  [8],  B<j  is  not 
assumed  to  equal  one,  but  is  derived  from  the  continuous 
model : 


fti  ♦  i 

Bd  =  <tB  ( ti  ♦  i  ,  T )  BdT  *  Btb  (1-$b  ) 


(3-60) 


where  B»l,  and  the  discrete,  augmented  control  input  matrix, 
fid  a  ,  is 

Bd •  -  [  Bd  0  0  0  0  ]T  (3-61) 

Because  the  cost  weighting  function  includes  a  cost 
weighting  term  on  the  pseudointegral  state,  Equation  (3-35) 
must  be  redefined  as 
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The  optimal  control  law  that  minimizes  the  cost  function  is 


u(ti  )  -  -£c*  (ti  ) Xa  (tt  )  +  Eyd 


-  [Gc  i  *  Gc  2  *  ] 


XB 

(ti ) 

XT 

p  (ti ) 

XT 

» (ti ) 

XT 

A  (tl  ) 

ti  ) 

+  Eyd 


(3-63) 


where  Gc  i 


a  [  Gcii*  Gcip*  Get**  Get  a*  ].  The  steady 
state  controller  gains  are  found  by  solving 


£e*  -  [u  +  &daT£ca£da  ]"1  [BdaTKcaOa  ] 

=*  Bd  [ki  i  ki  2  ki  3  ki  a  ki  a  ]  Oa  /  (U  +  Bd  2  ki  i  ) 


(3-64) 
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E  ■  [Gc  i  *  —  Gc  *  *  1 1  T /Kc  a  2  ]  Hi  *  +  Hi: 


where 


Ei  t 

Ei*  1 

•»  —  1 

0  : 

Bd 

*  o 

-I  : 

0 

ns 

Em 

£c 

-£r  i 

0 

(3-65) 


(3-66) 


Sea  is  the  steady  state  solution  to  the  backward  Riccati 
difference  equation,  and  the  solution  is  found  by  solving 
Equation  (3-50)  with  the  successive  values  of  JCca  equal  to 
one  another: 


Kc.  -  X.  +  -  [&».»&•£«  3TGc* 

where  the  elements  of  Kca  are  defined  as 


(3-67) 
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(3-68) 


Table  3-1  contains  the  PI  controller  gains  used  for  the 
sensitivity  and  robustness  analyses  (to  be  discussed  in 
Chapter  4) . 


TABLE  3-1 


Calculating  the  PI  Controller  Gains 


1.  Problem  Setup 


Tt  -  20 


tt  =  10 


.0001 


U  *  1 

Xn=  100 

Xa  a 

■  10 

2 .  Bd a  i 

4>a  ,  £ca  and  J1  Matrices 

Bda-  t  .97541 

o 

• 

o 

0.0 

0.0 

0.0  ] 

.95123 

0.0 

0.0 

0.0 

0.0 

0.0 

.99990 

.99990 

.48371 

0.0 

®e  * 

0.0 

0.0 

.99990 

.95158 

0.0 

0.0 

0.0 

0.0 

.90484 

0.0 

1.0 

-1.0 

0.0 

0.0 

1.0 

138.57 

-138.63 

-1.0092 

.21089 

37.367 

-138.63 

151.13 

62501 

624373 

-37.370 

Rea  c 

-1.0099 

62501 

6.25*10* 

6.25*10* 

-.00963 

.20405 

624373 

6.25*10* 

6.25*1010 

.69494 

37.367 

-37.370 

-.00946 

.69654 

47.099 

3.01*10- 

7  -10001 

-1.00*107 

-1.00*10® 

1.0000 

3.01*10“ 

7  -10001 

-1.00*107 

-1.00*10® 

1.8*10-® 

11  » 

0.0 

0.0 

-10001 

-1.00*10® 

0.0 

0.0 

0.0 

0.0 

-10.508 

0.0 

1.0252 

-500.02 

-5.00*10* 

-5 . 00*107 

.05000 

3.  PI  Controller  Gains 

£e  -  t  1.2422  -1.2921  -1.0252  -.49800 

E  -  1.074560 

.27436  ] 

Note:  c  is  the  distance  that  the  two  poles  are  moved  left 

of  the  origin  (see  Equation  3-39) 

is  calculated  from  Equations  (3-60) 

$•  is  calculated  from  Equations  (3-57) 

(3-59) 

Sea  is  calculated  from  Equation  (3-67) 
numerical  precision  that  is  required 
H  is  calculated  from  Equation  (3-66) 

£e  is  calculated  from  Equation  (3-64) 

E  is  calculated  from  Equation  (3-65) 


and  (3-61) 
through 

-  note  the 


I 

3.4  Reviewing  the  Cost  Weighting  Matrices 

This  section  will  re-evaluate  the  quadratic  cost 
function  in  both  the  proportional  gain  and  the  PI  trackers 
in  the  hope  of  selecting  a  set  of  weighting  matrices  that 
will  improve  the  robustness  of  the  trackers.  We  will  start 
with  the  more  complex  PI  tracker  cost  function,  and  reduce 
the  general  quadratic  cost  equation  to  the  present  form 
expressed  by  Equation  (3-45) .  After  reviewing  the 
relationship  between  the  elements  within  the  weighting 
matrices,  we  will  draw  the  analogues  applicable  to  the 
proportion  gain  tracker  cost  function. 

The  general  form  of  the  PI  cost  function  is 
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s  ( t*  ♦  i  ) 
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JCr  i  i  Krn 

x(  t«  ♦  l  ) 

g  ( tn  ♦  i ) 

Sr  l  a  T  Sr  a  a 

g(t«*i  ) 

To  generate  a  steady  state  control  law,  we  can  let  N 
approach  «•.  As  a  result,  the  matrix  will  not  effect  the 
overall  cost  function,  and  can  be  dropped  from  any  further 
consideration. 

The  £  cross  terms  have  two  origins.  First,  the  §_  term 
arises  from  the  coupling  between  the  control  input,  u(ti ) , 
the  pseudo-integral  vector,  g(ti )  ,  and  the  dynamics  state 
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vector,  x(ti ) .  This  coupling  arises  from  the  linear 
relationship  of  the  filter  model,  such  as, 

*k  ■  -  ( 1/ T )  Xk  +  Uk  (3-70) 

and  there  might  be  a  desire  to  put  a  quadratic  cost  on  the 
scalar  quantities  such  as  Xk  that  are  linear  combinations  of 
u,  q  and  &.  Second,  the  S  term  is  frequently  used  to  exert 
control  over  continuous-time  dynamics  for  the  entire  sample 
period  and  not  just  at  ti  [14] .  The  significance  of  this 
contribution  to  the  S.  term  depends  on  the  sample  rate.  As 
the  sample  rate  is  increased,  the  £  term’s  contribution 
approaches  zero.  For  the  controller  designs  in  this 
application,  it  is  assumed  that  the  sample  rate  will  be  high 
enough  to  make  the  S  terms  insignificant.  Without  the  S 
terms.  Equation  (3-69)  assumes  the  form  found  in  Equation 
(3-45)  . 

What  remains  to  be  evaluated  are  the  X*  and  U  cost 
weighting  matrices  required  for  the  PI  tracker  steady  state 
gain  functions.  The  &•  matrix  is  defined  as  a  positive 
semidefinite  weighting  matrix  that  assigns  a  cost  to  the 
tracker  error  between  xb  (ti )  and  xrr(ti).  Therefore,  there 
is  no  reason  to  weight  the  cross  product  terms  resulting 
from  the  target  velocity  or  acceleration,  and  their  elements 
are  set  to  zero  as  depicted  by  Equation  (3-62).  The  Xi  i 
elements  weight  the  present  position  error  between  the  beam 
and  target  and  Xss  weights  the  pseudointegration  value  of 


those  errors.  Non-zero  Xa a  allows  the  controller  to  have 
type  1  versus  type  0  properties,  which  are  desirable  in 
tracking.  If  Xa  a  is  smaller  than  Xu,  this  tells  the 
controller  to  place  a  greater  emphasis  on  minimizing  the 
integral  of  past  errors.  Therefore,  the  response  to  a  large 
positive  error  history  at  ti  would  be  a  large  negative 
control  input.  Unfortunately,  this  can  cause  a  very 
oscillatory  response  and  a  large  initial  overshoot  that 
could  drive  the  actuators  into  saturation  [8] .  If  this  were 
to  occur,  the  oscillatory  response  can  be  dampened  by 
increasing  Xg  a . 

The  U  matrix  places  a  penalty  on  the  control  energy 
expenditure,  and  U  must  be  positive  definite.  Because  we 
are  limiting  the  design  to  one  dimension,  U  becomes  a 
positive  scalar.  This  leaves  only  three  cost  weighting 
elements,  Xi i ,  Xa a  and  U,  to  be  evaluated.  If  one  of  the 
elements  is  set  to  a  constant  value,  then  a  robustness 
analysis  can  be  conducted  over  the  range  of  the  remaining 
two  elements.  In  other  words,  the  selection  of  alternate 
weighting  matrices  that  have  the  best  robustness 
characteristics  does  not  require  us  to  evaluate  the  the 
three  by  three  cost  weighting  matrix  depicted  by  Equation 
(3-69) ,  but  instead  requires  us  to  confine  our  attention  to 
a  single  relationship  with  only  two  unknown  variables. 

These  findings  can  be  extended  to  the  proportional  gain 
tracker.  Since  the  proportional  gain  tracker  does  not  have 
a  pseudointegrator,  it  does  not  have  a  pseudointegral 


weighting  element.  Therefore,  if  we  set  U  to  a  constant, 
then  we  only  need  to  vary  Xi 1  to  determine  the  best  set  of 
cost  weighting  matrices. 

The  method  of  selecting  a  set  of  weighting  matrices  is 
based  on  the  simple  sum  of  the  quadratics  on  all  the 
quantities  of  interest  at  each  sample  period.  Therefore, 
any  deviations  from  zero  are  penalized.  If  we  had  to  design 
a  system  with  specified  performance  limitations,  we  would 
pick  the  initial  cost  weighting  elements  as  one  over  the 
square  of  the  maximum  allowable  value  [14].  From  there,  the 
controller's  cost  weighting  matrix  would  be  tuned  according 
to  the  results  of  a  performance  analysis.  Without 
prescribed  design  specifications,  the  cost  weighting 
elements  will  be  selected  to  minimize  the  RMS  tracker  error 
which  will  be  calculated  during  the  initial  performance 
analysis  (refer  to  Section  4.4). 

There  are  several  alternative  methods  of  selecting  cost 
weighting  matrices.  First,  with  a  controller  feedback 
sample  period  of  one  second,  it  might  be  desirable  to 
include  the  cross  terms  in  order  to  provide  better  control 
over  the  entire  sample  period.  Second,  the  weighting 
matrices  can  be  selected  on  the  basis  of  their  impact  on  the 
closed-loop  system's  pole  placement.  The  effects  this  can 
have  on  stability  and  loop  shaping  can  be  observed  by 
evaluating  the  effects  different  controller  gains  could  have 
on  the  controller  transfer  function  (see  Equation  (6-13)). 
Two  formal  techniques  which  exploit  this  are  implicit  model 


following  [2,16,19]  and  the  LQG/LTR  Dual  Method  of 
Kwakernaak  and  Sivan  [10] .  Further  discussions  on  these 
techniques  has  been  included  in 
Appendix  A. 

3.5  The  Multiple  Model  Adaptive  Controller 

The  previous  two  controller  designs  had  difficulty 
tracking  a  target  transitioning  from  a  benign  trajectory  to 
an  evasive  flight  path  because  of  the  dynamics  noise 
strength  parameter  uncertainty  induced  by  the  increased 
maneuvering.  That  is,  if  a  non-adaptive  Kalman  filter  were 
tuned  to  a  target  flying  straight  and  level,  it  will  have 
great  difficulty  tracking  the  same  target  as  the  target 
tries  to  evade  the  tracker  with  a  series  of  high-g  jinking 
maneuvers.  The  reason  is  that  the  filter's  parameters  such 
as  the  target's  acceleration  time  constant  and  the  dynamics 
noise  strength  will  not  be  representative  of  the  new 
target's  trajectory.  One  method  to  improve  a  system's 
performance  is  to  construct  a  multiple  model  adaptive 
controller  (MMAC) .  This  section  will  start  out  by 
developing  the  general  MMAC  model  and  then  tailor  it  to  the 
tracker  model  that  is  used  throughout  the  thesis.  The 
uncertain  parameters  of  initial  interest  are  those 
associated  with  the  target.  Thus,  the  multiple  model 
structures  entail  replications  of  the  Kalman  filter  or  of 
the  entire  controller. 

3.5.1  Derivation  of  the  Multiple  Model  Adaptive 


The  MMAC  is  based  on  the  premise  that  the 


robustness  of  a  controller  can  be  improved  if  the  controller 
can  provide  on-line  estimation  of  the  uncertain  parameters 
defined  is  a ( t) .  The  uncertain  parameters  can  affect  any 
or  all  of  $,  gd ,  H,  Qd  and  R,  and  thereby  induce  erroneous 
estimates  of  the  state  and  the  error  covariance  in  a  non- 
adaptive  filter.  As  a  result,  we  may  not  be  able  to  design 
a  single  controller  model  exactly.  Instead,  we  may  find  it 
possible  to  design  an  adaptive  estimator /controller . 

One  means  of  generating  an  adaptive  stochastic  feedback 
controller  is  by  designing  a  MMAE-based  adaptive  controller 
with  adaptive  controller  gains,  Gc* [ti ,g(ti ♦ ) ]  (see  Figure 
3-5a) .  The  concept  is  that  the  MMAE  provides  a  state 
estimate,  £(ti+),  and  an  uncertain  parameter  estimate, 
S(tiM,  based  on  the  past  measurement  history,  Zi  .  The 
adaptive  controller  gain  is  calculated  from  the  uncertain 
parameter  estimate  (gains  are  precomputed  as  a  function  of 
assumed  parameter  values,  and  then  those  functions  are 
evaluated  based  on  SHti*)),  and  the  control  input,  u(tt ) ,  is 
calculated  from  the  product  of  the  state  estimate  and  the 
adaptive  controller  gain,  as: 

H(ti)  -  -Gc*  t  ti  ,I(tt  )J*(ti*  )  (3-71) 

If  the  uncertain  parameters  were  confined  to  the 
dynamics  driving  noise  and  the  measurement  corruption  noise 
statistics,  then  GpMti)  is  not  dependent  on  S(tt ) ,  and  a 
non-adaptive  controller  gain,  GcoMti),  is  appropriate  (see 


Figure  3-5.  Adaptive  Stochastic  Feedback  Controller: 

<a)  With  an  adaptive  feedback  gain.  (b)  With  non-adaptive 
gain. 


Figure  3-5b) .  On  might  also  choose  to  implement  a  nominal 
gain  Geo (ti )  rather  than  employ  an  adaptive  gain,  using 
adaptivity  only  to  improve  the  accuracy  of  the  state 
estimate.  Last,  one  might  assume  steady-state  conditions 
and  use  a  steady-state  controller  gain,  such  as  Gc*(t)  or 
Geo*. 

By  extending  the  concept  of  an  adaptive  stochastic 
feedback  controller,  we  can  generate  the  multiple  model 
adaptive  controller  (see  Figure  3-6)  [1,4,13,14].  This  is 

done  by  replacing  the  bank  of  Kalman  filters,  used  to 
estimate  the  state  and  uncertain  parameters,  with  a  bank  of 
LQG  controllers.  Each  controller  is  based  on  a  particular 
set  of  parameter  values,  ait  (ti  )  ,  where  k  identifies  the 
controller  within  the  bank.  Within  each  controller,  a  state 
estimate,  kn(tiM,  is  computed  by  a  Kalman  filter  based  on 
the  parameter  value,  ait  ,  and  is  multiplied  by  the 
appropriate  steady-state  gain,  Gc*(an),  based  on  the  same 
parameter  assumption.  The  result  is  an  optimal  control  law 
conditioned  on  &k  being  the  true  parameter: 


liu(ti)  -  -fie*  (a>)gic  (ti  ♦  )  (3-72) 


The  adaptive  control  is  generated  by  adding  the 
probabilistically  weighted  Uk  (ti  )  values;  as  shown  at  the 
right  most  summing  junction  in  Figure  3-6. 
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Multiple  Model  Adaptive  Controller 
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One  way  to  develop  the  HMAC  mathematically  is  first  to 
develop  the  multiple  model  adaptive  estimator  using  a 
Bayesian  approach,  and  to  draw  the  mathematical  analogy 
between  the  adaptive  feedback  controller  based  on  such  an 
estimator  and  the  MMAC .  The  Bayesian  state  estimate  can  be 
expressed  as  the  expected  value  of  the  state  conditioned  on 
the  measurement  history  Zi : 


(3-73) 


*(tiM  *  E{&(ti  )  I  Zi  }  -  1-fx  <  t  >  i  2  (ll£i  )dl 


where  i  is  a  dummy  variable  of  integration  of  x(t)  ,  and 
fs  <  t  >  i  2  (£.l£i  )  is  the  conditional  probability  density 
function  of  x(t) .  If  we  let  a(t)  denote  the  vector  of 
uncertain  parameters,  and  assume  &(t)  can  take  on  any  value 
in  the  continuous  range  of  A,  Equation  (3-73)  can  be 
rewritten  to  include  a(t) : 


k(ti*)  *  £  f*  <  t  > ,  •  i  z  (£.,al2ii  )da  di. 


(3-74) 


where  g  is  the  dummy  variable  used  to  integrate  &( t)  over 
the  range  of  A.  According  to  Bayes'  rule,  the  conditional 
probability  density  function  can  be  expressed  as 


fx(t,,«.z  <i,slZi  )  “  f*<t)i«,z(.£.la,Zi)f«!2(a!2i.i) 


(3-75) 


v.  / 

•* 


65 


if 


I 


6C 


$ 


The  first  density  function  on  the  right  hand  side  of 
Equation  (3-75)  is  Gaussian,  with  a  mean  of  5UtiM  and  a 
covariance  of  P(tiM  as  computed  by  a  Kalman  filter  based 
upon  a  particular  a.  The  second  term  is  the  conditional 
probability  of  a  conditioned  on  the  measurement  history. 
According  to  Bayes*  rule,  it  can  be  expressed  as 

f  a  t  2  (<2.1  Zt  )  -  faiz(i), 2(1-1)  ( Q  I  £.1  ,  Z  i  -  i  ) 


is 


f  a  ,  z  (  1  )  I  2  (  1  -  l  )  (a  ,  I  Z l  -  l  ) 

f  z  <  1  )  I  2  (  i  -  i  >  (£i  I  Z.1  -  1  ) 

f  z  (  i  >  i  a  .  2  (  i  -  i  >  (£i  I  a,  Zi  -  i  )  f  a  i  2  (a  !  Zi  ) 

f z  <  i  >  i  a  ,  2  <  i  -  i  >  (ii  la,  Zt  -  i  )  fa  i  z  (a  I  Zt  )  da 

(3-76) 
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where  f  z  <  i  >  i  a  .  z  <  i  -  i  >  (£i  I  a,  £i  -  i  )  is  Gaussian ,  with  a  mean  of 
g(ti)£(ti")  and  a  covariance  of  [H(ti  )  P(ti  -  )||T  (ti  )  +  R(ti)]. 
Equation  (3-76)  could  conceptually  be  evaluated  recursively 
starting  from  fa  (q) .  By  combining  Equations  (3-74)  with 
(3-75)  and  switching  the  order  of  integration,  we  can 
generate  the  state  estimate  from 


*(ti  ♦  ) 


Ja  >- 


I  fi<t>ia.z(£.lg,£i)d£  fa:z(glZt)da 


(3-77) 


Unfortunately,  the  double  integration  will  make  this 
approach  computationally  impractical,  so  we  let  the 
uncertain  parameter  vector  assume  only  a  finite  number  of 
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values,  such  as  the  discrete  vector  set  {ai ,  a 2 ,  ...  a* l , 
where  K  is  the  number  of  elemental  filters  in  the  MMAE ,  and 
eventually  the  number  of  controllers  within  the  MMAC .  Thus, 
the  integration  is  replaced  by  a  finite  summation.  This 
becomes  the  foundation  of  the  multiple  model  structure, 
where  each  controller  within  the  bank  of  controllers  is 
based  on  an  actual  discrete  parameter  value. 

To  adopt  the  Bayesian  approach  to  a  multiple  model 
structure,  an  a  priori  density  function  is  required  for 
3*  ( to  )  : 

K 

fa  <  t  «  )  (a)  =  I  pic  ( to  )  5  (a  -  ak  )  (3-78) 

It  *  1 


where  pk(to)  is  the  probability  a  assumes  the  value  an  at 
time  to.  Therefore,  the  hypothesis  conditional  probability, 
pn (ti ) ,  is  a  recursive  relationship  expressed  as 


Pk  (ti  )  -  prob(a=a.k  IZ(ti  )  *Zi  } 

f  z  (  1  >  1  a  ,  2(1-1)  ( £.1  I  ak  ,  Zl  -  1  )  Pk  ( 1 1  -  1  ) 


(3-79) 


I  fz  <  1  )  !  a  ,  Z  <  1  -  1  )  (£.1  I  Su  »Z.I  -  1  )  PJ  ( tl  -  1  ) 
J  - 1 


where  f z  <  1  >  i  a  ,  2  <  1  -  1  >  (£,i  I  &k  ,  £1  - 1  )  can  be  evaluated  as 


f  z  (  1  >  1  a  .  z  <  1  -  1  )  (£.t  I  a>  ,  Zi  -  1  ) 


(2n)-/>  I Ak  (ti  ) 1 1 / » 


exp ( • I 


(•I  -  l-Hr>T  (ti  )^k  (ti  )-‘  &  (ti  ) 


(3-80) 
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and 


I>  (ti  )  is  the  residual;  j>  (ti  )  -  2.1  -  gic  (ti  )gk  (ti  -  ) 
&k(ti)  »  fik  (ti  )£).  (ti- )Hk  (ti  )  +  Rk  (ti  ) 


m  is  the  number  of  measurements  at  time  tt ;  dimension 
of  2)1  (ti  ) 


Therefore,  the  conditional  mean  is  the  sum  of  the 


probabilistically  weighted  state  estimates  from  each 


controller  within  the  bank: 


*(tt*  )  *  E(x(ti  )  IZ(ti  )»£t  \ 


I  *k  (ti  ♦  )pk  (ti  ) 
k  -  1 


(3-81) 


Note  how  equations  (3-79)  and  (3-81)  are  closely  related  to 


equations  (3-76)  and  (3-77)  respectively. 


The  MMAC  feedback  control  can  be  analogously  derived, 


and  is  the  probabilistically  weighted  sum  of  jjk  (ti  )  as 


expressed  by 


u(tt  ) 


I  JAk  ( ti  )  Pk  ( 1 1  ) 
k  «  l 


I  — £e  *  (a,k  )£k  (ti  ♦  )pk  (ti  ) 
k  >  1 


(3-82) 


3.5.2  The  Multiple  Model  Adaptive  Tracker  Since  this 
is  the  first  implementation  of  a  MMAC  for  this  application 
and  is  to  be  viewed  as  a  feasibility  demonstration,  the 
design  will  be  simple.  The  MMAC  will  be  based  upon  only  one 
uncertain  parameter,  discretized  into  three  possible  values, 
and  thus  it  will  be  composed  of  a  bank  of  three  controllers. 
The  purpose  of  the  MMAC  is  to  allow  the  tracker  to  adapt  to 
changes  in  the  best  modeled  value  of  target  dynamics  driving 
noise  strength,  Qt ,  because  the  tracker  has  no  control  over 
Qt ,  and  since  mis-modeling  of  Qt  leads  to  large  RMS  errors. 
Each  controller  will  be  tuned  to  a  different,  quantized 
value  of  Qt k .  For  this  application,  each  controller  will  be 
tuned  to  one  of  the  following  uncertain  parameter  values: 

Qti  *  0.01  Qt  i  s  0  •  1  Qt  i  =  1 . 0 

These  three  values  approximately  cover  (for  modeling 
purposes  only)  the  range  of  trajectories  from  the  straight 
and  level  flight,  to  the  maximum  manned-vehicle  g-limit  of 
lOg's,  assuming  a  specific  range  to  target  or  a  defined 
relationship  between  real  maneuvers  and  maneuvers  as  seen  in 
the  detector  array  plane. 

Because  Qt  is  selected  as  the  uncertain  parameter  and 
it  is  not  used  to  calculate  the  steady-state  controller 
gain,  we  can  use  the  non-adaptive  steady-state  controller 
gain,  Qeo*.  If  we  had  to  calculate  an  adaptive  controller 
gain,  such  as  would  be  the  case  if  the  target  acceleration 

* 
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tine  constant,  tt ,  were  the  uncertain  parameter,  then  each 
controller  gain,  Gjc  *  [aic  { ti  )  ]  ,  would  be  calculated  from  the 
same  steady-state  controller  gain  equations.  Equations 
(3-57)  through  (3-67) ,  but  each  controller  gain  function 
would  be  tuned  to  its  specific  uncertain  parameter  value. 

Each  controller  within  the  MMAC  structure  operates 
independently  of  one  another,  and  the  embedded  k-th 
elemental  Kalman  filter  calculates  the  residual,  ru (ti ) ,  the 
covariance  of  the  measurement,  Ak(tic),  and  the  best  estimate 
of  Xk (ti  ) .  The  residuals  and  covariances  of  the  measurement 
from  all  of  the  elemental  filters  are  used  to  calculate  the 
probability,  pk (ti ) ,  which  is  used  to  weight  each  controller 
(see  Equations  (3-79)  and  (3-80)).  Then,  the 
probabilistically  weighted  control  inputs,  uk  (ti )pk (ti ) ,  are 
summed  according  to  Equation  (3-82) .  The  controllers  start 
out  with  equal  probabilistic  weighting.  That  is,  each 
controller  is  initialized  with  an  equal  probability,  pk(to) 
■  1/K;  the  constant  K  represents  the  number  of  controllers 
within  the  MMAC  structure  (for  this  application,  K  =  3  ) . 

3 . 6  Summary 

This  chapter  developed  the  target  model  and  all  the 
controllers  to  be  used  directly  or  in  support  of  this 
thesis.  The  target  was  modeled  through  a  first  order  Gauss- 
Markov  acceleration  process,  where  the  state  estimate 
(involving  position,  velocity  and  acceleration  variables) 
and  covariance  error  are  provided  by  a  Kalman  filter.  The 
first  controller  design  developed  the  proportional  gain 


regulator  and  tracker  equations.  But  because  the 
proportional  gain  controller  does  not  have  type-one 
characteristics  that  allow  rejection  of  unknown,  constant 
disturbances  that  arise  from  linearized  models,  the  PI 
controllers  were  developed.  Still,  the  performance  of  the 
tracker  was  sub-optimal.  The  next  sections  looked  at 
improving  the  controllers '  robustness  by  evaluating 
alternative  weighting  matrices  and  using  the  PI  tracker  in  a 
multiple  model  adaptive  controller.  In  all,  the  final 
objective  is  to  design  the  optimal  controller  which  can 
place  the  particle  beam  on  a  maneuvering  target. 
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IV.  THE  ANALYSES 

Before  a  design  can  be  considered  complete,  it  must  be 
developed  and  its  performance  evaluated  against  some 
prespecified  standard  or  baseline.  This  chapter  discusses 
the  different  performance  analyses  used  to  evaluate  the  PI 
controller  and  also  to  develop  and  evaluate  the  adaptive  PI 
controller.  Both  controller  designs  use  the  Meer  filter  to 
estimate  the  beam  (controlled)  state,  and  the  Kalman  filter 
to  estimate  the  target  (reference)  state.  The  results  of 
the  performance  analyses  are  listed  in  Chapter  5.  The  first 
section  explains  why  Monte  Carlo  simulation  was  selected 
over  covariance  analysis.  The  next  section  addresses  the 
software  packages  SOFE  and  SOFEPL  which  were  used  to 
generate  the  Monte  Carlo  simulation  and  error  statistics. 

The  following  section  provides  an  in-depth  review  of  the 
parameters  used  in  the  Meer  filter  for  beam  location 
estimation  and  the  Kalman  filter  for  target  state 
estimation.  The  last  section  discusses  the  actual  analyses 
that  are  used. 

4.1  The  Method  -  Monte  Carlo  Simulation 

The  particle  beam  estimator/controller  problem  is 
modeled  as  a  stochastic  process  and  the  performance  can  be 
best  evaluated  by  observing  the  statistical  behavior  of  the 
error  processes.  Specifically,  we  are  interested  in 
evaluating  the  controller's  tracking  error  statistics,  and 
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it  is  assumed  that  the  embedded  estimators  are  already  tuned 
for  the  appropriate  dynamics  driving  noise  strength,  Q,  and 
measurement  corrupting  noise  covariance  (spread  dispersion 
of  the  beam  in  the  Meer  filter),  R.  Regardless,  it  will 
also  be  necessary  to  evaluate  the  Kalman  and  Meer  filter's 
state  estimation  error  statistics,  to  support  some  of  the 
findings  from  the  sensitivity  and  robustness  analyses.  The 
performance  of  a  controller  can  often  be  analyzed  by  either 
covariance  analysis  or  Monte  Carlo  simulation.  The  more 
efficient  method  is  the  covariance  analysis,  which  requires 
only  one  software  run  to  generate  the  time  history  of  the 
estimation  error  covariance,  P«  (ti ) ,  or  other  pertinent 
statistics.  Unfortunately,  covariance  analysis  cannot  be 
applied  to  the  multiple  model  adaptive  controller  because  a 
proper  covariance  analysis  is  limited  to  linear  stochastic 
controllers  that  use  prespecified  measurement  update  times. 
The  adaptive  mechanisms  in  both  the  controller  and  in  the 
Meer  filter  violate  the  linearity,  and  the  elemental  Snyder- 
Fishman  filters  have  a  varying  sample  rate  that  is  not 
prespecif iable . 

Therefore,  the  less  efficient  Monte  Carlo  simulation  is 
used  to  evaluate  the  particle  beam  estimator/controller. 

The  Monte  Carlo  simulation  is  a  complete  computer  simulation 
that  requires  numerous  simulation  runs  to  generate  enough 
samples  of  the  error  process  to  approximate  the  filter's 
true  statistics  with  sample  statistics  [12] .  For  this 
application,  Zicker  [17,27]  found  that  the  sample  statistics 


sufficiently  converged  by  200  runs,  with  each  run  lasting 
100  seconds.  Zicker  also  found  that  the  system  transients 
dissipated  by  50  seconds  into  the  simulation.  Therefore, 
these  numbers  have  become  the  baseline  for  future  Monte 
Carlo  simulations  and  the  statistics  of  the  error  processes 
are  averaged  over  the  last  50  seconds.  Because  the  last  50 
seconds  should  represent  a  steady-state  condition,  the  time- 
average  of  the  error  process  statistics,  such  as  the  average 
root-mean-square  tracker  error,  can  be  used  as  a  compact 
index  of  performance,  and  can  be  used  to  identify  trends 
within  any  of  the  analyses. 

The  Meer  or  Kalman  filter's  state  estimation  error 
statistics,  which  are  used  to  evaluate  the  filter's 
operation,  are  the  mean  and  the  standard  deviation  of  the 
error  between  the  truth  and  filter  states.  They  will  be 
calculated  and  plotted  for  each  controller  feedback  sample 
period,  ti ,  between  50  and  100  seconds.  For  example,  the 
derivation  of  the  Meer  filter's  beam  position  state  estimate 
error  statistics  are  as  follows:  the  error  between  the 
true  beam  position  and  the  Meer  filter  estimate  is 

ea  (ti  ,n)  ■  xta  (ti  ,n)  -  fta  (ti  ,n)  (4-1) 

where  50  <.  ti  <.  100  and  n  is  the  run  number.  Thus,  the 
mean  error  is 

N 

«• (ti )  -  [l/N] -I  e>  (ti  , n) 


(4-2) 


where  N  is  the  total  number  of  runs  (N»200) .  The  variance 
and  the  standard  deviation  of  the  error  are  [5,7] 

M 

v.  (ti  )  «  [1/  (N-l)  ]  •  I  {e»Mti,n)|  -  [N/ (N-l)  ]  •  @b  2  ( ti  )  (4-3) 

■  =«  1 

Ob  (ti  )  -  (1  +  1/[4(N-1)]  }  •  (Vi  (ti  )  (4-4) 

Equation  (4-3)  calculates  the  unbiased  estimate  of  the 
variance  and  Equation  (4-4)  is  an  approximation  that 
produces  an  unbiased  estimate  of  the  standard  deviation.  A 
full  derivation  of  the  exact  equation  and  the  limitations 
that  apply  to  the  approximation  can  be  found  in  reference 
[6] .  Because  these  statistics  are  calculated  for  the 
steady-state  condition,  it  is  desirable  for  §b  (ti  )  to 
approximate  zero  (i.e.,  it  is  desirable  for  the  bias  to  be 
negligibly  small) ,  and  for  oi (ti )  to  approach  a  steady  state 
value  from  ti »50  to  100  seconds.  Deviations  from  these  two 
conditions  will  indicate  the  filter  is  operating  improperly. 
Note  how  the  statistics  are  calculated  for  a  finite  sample 
population,  and  only  when  we  have  a  large  enough  sample  can 
we  be  sure  the  sample  statistics  closely  approximate  the 
true  statistics. 

The  root-mean-square  (RMS)  tracking  error,  RMS. (ti ) , 
and  the  time-averaged  RMS  tracker  error,  RMS. ,  are  the 
statistics  used  to  evaluate  the  performance  of  the 
controller  during  the  different  analyses.  The  error  for  the 


tracker  problem  is  defined  as  the  difference  between  the 
true  target  state  and  the  true  beam  position  state: 

et  <ti  ,n)  »  xt«  (ti  ,n)  -  xtxp  (ti  ,n)  (4-5) 

The  errors  from  each  run  are  averaged  and  a  mean  error, 
m*  (ti ) ,  and  a  standard  deviation  of  error,  o*  (ti ) ,  are 
calculated  for  each  ti : 


N 

m*  (ti  )  =  [1/N]  *  Z{xt«(ti,n)  -  xtTp(ti,n)I  (4-6) 

a  *  1 

o.  (tt  )  -  11  +  1/  [4  (N-l)  ]  l 

M 

•  I  [1/  (N-l)  ]  •  Z  f  e» 2  (ti  ,n)  }  -  [N/  (N-l  )]*nu2(ti)|1/2 

«*i 

(4-7) 

The  root-mean-square  (RMS)  was  selected  as  the 
statistical  performance  parameter  to  be  used  to  evaluate  the 
controller.  The  RMS  error  can  be  generated  for  each  sample 
period  by  using  the  equation: 

RMS*  ( tt  )  -  [m*Mti  )  +  o. 2  (ti  )]»/*  (4-8) 


n 


The  time-averaged  RMS  error  is  calculated  by  averaging  each 
RMS  error,  RMS*  (ti ) ,  from  50  to  100  seconds: 


too 


RMS.  -  [1/51]  *Z  RMS.(ti) 


(4-9) 
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Once  again,  because  the  RMS  error  statistic  is  calculated 
for  steady-state  condition,  RMS*  (ti )  should  appear  as  a 
constant  value  with  only  minor  deviations.  A  controller 
with  an  increasing,  ramp-like  RMS  error  would  indicate  that 
the  filter  is  unstable.  Anytime  an  undesirable  RMS  error  is 
found,  the  RMS  error  will  be  analyzed  to  see  if  the  problem 
is  primarily  due  to  bias-like  characteristic  (i.e.,  an 
increasing  mean)  or  fluctuations  (i.e.,  an  increasing  o* ) . 

4.2  The  Tools  -  SOPE  and  SOFEPL 

The  performance  analyses  are  generated  from  tailored 
software  packages  called  SOFE  and  SOFEPL  [7,21].  SOFE, 
which  is  short  for  Simulation  for  Optimal  Filter  Evaluation, 
is  a  general  purpose  Monte  Carlo  simulation  program  designed 
for  evaluating  systems  that  use  Kalman  filters.  It  contains 
both  a  truth  model  and  a  filter  to  be  evaluated.  The  truth 
model  is  described  by  a  set  of  stochastic  differential 
equations  that  emulate  the  "real-world”  system  dynamics. 

The  filter  is  also  described  by  a  set  of  differential 
equations  representing  the  system's  state  propagation  within 
a  Kalman  filter,  and  updates  due  to  measurements.  SOFE 
allows  the  user  to  specify  the  measurement  format  of  the 
filter  update,  and  if  desired,  to  include  feedback  control. 
SOFEPL,  which  is  short  for  SOFE  Plotter,  is  a  post¬ 
processing  routine  that  can  be  programmed  to  perform  various 
statistical  functions,  such  as  calculate  the  RMS  error  time 
history,  and  then  plot  the  results.  Both  SOFE  and  SOFEPL 


were  developed  by  Stanton  H.  Musick  and  others  of  the  Air 
Force  Avionics  Laboratory. 


SOFE  had  to  be  modified  internally  before  it  could  be 
used,  because  the  original  version  of  SOFE  was  limited  to 
Kalman  or  extended  Kalman  filter  applications.  In  doing  so, 
Meer  [18]  changed  SOFE's  code  to  accommodate  the  Snyder- 
Fishman  estimators  and  the  MMAE  structure  for  indicating  the 
likelihood  of  each  photoelectric  event  being  due  to  signal 
versus  noise.  The  problem  of  adapting  SOFE  to  this 
application  did  not  affect  the  propagation  or  update 
equations,  but  lay  in  generating  Poisson  distributed  random 
signal  and  noise  events,  and  then  using  these  events  as 
measurement  update  for  the  Meer  filter.  The  standard  Kalman 
filter,  such  as  found  in  the  original  version  of  SOFE, 
updates  the  filter  at  regular  discrete  time  intervals, 
whereas  the  Meer  filter  updates  only  when  it  receives  a 
signal  event  [27]  . 

Because  SOFE  was  intended  for  general  applications,  it 
allows  the  user  to  tailor  SOFE  to  the  user's  needs  through 
nine  user-defined  (i.e.  user-written)  subroutines.  These 
subroutines  allow  the  user  to  specify:  the  truth  model’s 
differential  equations,  the  filter  model's  differential 
equation,  truth  model  disturbances  (optional) ,  the 
projection  matrix  of  the  state  vector  onto  the  measurement 
vector  ( U.)  ,  the  measurement  corruption  noise  covariance 
matrix  (&) ,  the  system  dynamics  matrix  (often  composed  of 
partial  derivatives)  (F) ,  dynamics  driving  noise  strength 
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matrix,  (Q) ,  the  initial  filter  covariance  matrix  (Po ) ,  the 
initial  conditions  for  each  run,  a  prespecified  trajectory 
(optional) ,  impulsive  control  (optional) ,  and  other 
calculations  required  during  the  measurement  update  such  as 
feedback  control.  Meer's  [18]  alterations  to  the  code 
incorporated  several  more  user-defined  subroutines  that  were 
used  to  implement  the  Meer  filter.  Zicker  [27]  followed 
Meer  and  added  the  "Merge"  method  of  filter  pruning  and  the 
first  controller  design.  Moose  and  Jamerson  [8,20]  designed 
the  follow-on  controllers  and  Jamerson  implemented  the 
Gauss-Markov  acceleration  model  for  the  target.  Zicker, 
Moose  and  Jamerson 's  changes  are  limited  to  the  user-defined 
subroutines  and  do  not  affect  the  executive  routine  of  SOFE. 
The  SOFE  source  code  that  Meer  altered  and  that  is  used 
throughout  the  research  is  the  FORTRAN  4,  May  1982  version 
of  SOFE. 

The  conceptual  operation  of  the  modified  SOFE  program 
is  shown  at  a  macro  level  in  Figure  4-1.  The  fixed-sample- 
period  controller  requires  that  the  discrete  feedback 
control  and  the  Meer  filter  propagation  cycle  operate  at  the 
same  discrete  sample  rate,  because  the  feedback  law  requires 
the  most  current  particle  beam  state  estimate.  As  mentioned 
before,  the  Meer  filter's  measurement  update  occurs  whenever 
an  event  is  observed.  For  simplicity,  the  time  to  the  next 
signal  and  noise  events  are  pre-calculated  from  the  signal 
and  noise  rate  parameters,  and  the  integrator  is  told  in 
advance  where  to  stop  for  the  Meer  filter  measurement 
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update.  A  diagram  of  the  top-down  structure  of  the  modified 
SOFE  is  presented  in  Appendix  B. 

The  propagation  of  the  truth  states,  the  filter  states 
and  the  filter-computed  error  covariances  is  accomplished  by 
integrating  the  differential  equations  from  the  current  time 
forward  to  some  specified  time.  This  time  can  be  the 
integration  step  size,  the  time  for  the  next  Meer  filter 
update  or  the  time  for  the  next  Kalman  filter  update.  The 
state  differential  equations  have  the  general  form 

*(t)  »  £[s<t) ,u(t) ,tj  +  G(t)w(t)  (4-10) 

where  w(t)  is  a  zero-mean  white  Gaussian  noise  of  strength 
fit-  SOFE  solves  Equation  (4-10)  by  first  using  a  fifth  order 
Kutta-Merson  integrator  to  solve  the  deterministic 
differential  equation 

*(t)  -  tts(t) ,n(t) ,t]  (4-11) 

The  stochastic  term  is  added  upon  the  completion  of  each 
integration  through  the  expression 

»  -  S>  +  GAUSS (O.Od )  (4-12) 

where  xd  is  the  deterministic  solution  to  Equation  (4-11) , 

25 s  is  the  stochastic  solution,  and  the  expression 

GAUSS (O.Od )  is  a  randomly  generated  vector  term  of  zero  mean 


and  a  covariance  of  Q,d .  The  dynamics  driving  noise 
strength,  Qd  (the  second  moment  of  the  equivalent  discrete 
time  noise  representing  G(t)w(t)),  is  calculated  from 

fti*i 

Qd  -  $(ti»i  ,t)G(t)Q(t)Qt  (t)$*  (ti  ♦!  ,T)dT  (4-13) 

Jti 

For  both  the  Kalman  and  Meer  filters,  Qd  is  scalar,  and  the 
integral  in  Equation  (4-13)  is  solved  as 

Qd  ■  GaQr [l-exp(— 2at/T) ] /2  (4-14) 

where  t  is  the  target  or  beam  time  constant,  G  **  1 ,  and  at 
is  the  sample  period  over  which  the  the  noise  is  injected 
into  the  Equation  (4-12) . 

The  measurement  noise  is  incorporated  into  the 
measurement  update  in  a  similar  manner.  The  general 
measurement  equation  for  either  the  Kalman  filter  or  for  the 
elemental  Snyder-Fishman  filter  (which  receives  a 
hypothetical  signal-induced  event)  is 

l(ti  )  »  &[£»  ( t*  ) , ti  ]  +v(ti)  (4-15) 

where  £(ti  )  is  the  sum  of  the  measurement  function  of  the 
true  states,  ( ti  ) , ti  ] ,  and  the  measurement  corruption 

noise,  y(tt ) .  SOFE  injects  the  measurement  noise  by  the 
following  relationship: 
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where  R  is  the  covariance  of  the  target  sensor  measurement 
noise  v(ti ) ,  or  the  beam  dispersion  for  the  beam  signal- 
induced  photodetector  event.  If  the  event  is  signal- 
induced,  Equation  (4-16)  gives  the  actual  location  of  the 
event,  not  a  "noise-corrupted"  measurement  of  it.  If  the 
event  in  noise-induced.  Equation  (4-16)  does  not  apply  and 
the  location  of  the  noise-induced  event  is  simulated  by 
using  a  uniformly  distributed  mapping  function. 

SOFE  updates  the  Kalman  filter  at  regular  time 
intervals  since  the  Kalman  filter  has  a  fixed  sample  period. 
This  is  not  true  of  the  Meer  filter  because  the  Meer  filter 
update  is  based  on  the  arrival  of  the  next  signal  or  noise 
event.  Thus,  Meer  [18]  had  to  modify  the  SOFE  source  code 
so  that  SOFE  would  simulate  a  varying  sample  period  (i.e., 
the  time  to  the  next  signal-  or  noise-induced  event)  as  a 
Poisson  time  process.  The  equation  which  will  calculate  the 
the  varying  sample  period  is  derived  from  the  Poisson 
process  density  function, 

P(y)  =  -{ [X (t) t]» /y ! I *exp [-X (t) t]  (4-17) 

where  p(y)  is  the  probability  that  y  events  occurred  within 
time  period  t,  and  £(t)  is  the  mean  signal  arrival  rate. 
Because  we  are  interested  in  generating  a  random  sample 
period  with  the  next  event  occurring  at  the  completion  of 


the  sample  period,  y  is  set  to  zero  (i.e.,  no  events  occur 
during  the  sample  period),  and  the  random  sample  period,  t, 
is  calculated  from  the  inverse  mapping  function  of  Equation 
(4-17)  [27] 


t  -  -i(t) -ln[p(0)] 


(4-18) 


where  p(0)  is  defined  as  the  probability  that  t  amount  of 
time  has  passed  before  an  event  arrives.  This  probability 
has  a  range  of  0£p(0)£l.  Because  we  are  interested  in 
generating  samples  of  ti ,  p(0)  is  selected  randomly  from  the 
range  of  p(0),  and  this  process  is  used  to  simulate  the 
arrival  time  of  the  next  signal-  or  noise-induced  event. 

The  mean  signal  arrival  rate  is  defined  as 


is  (t)  -  n/(tr-to  ) 


(4-19) 


where  n  is  the  number  of  signal-induced  events  expected  over 
the  duration  of  the  simulation  run,  (tF-to).  The 
relationship  between  the  mean  signal  arrival  rate  and  the 
expected  signal  rate  parameter.  &  (t ,£,&(t) ) ,  can  be  shown 
by  integrating  Equation  (2-1) , 


ee 

is  (t)  »  is  (ti  ,g)da 

— e# 


A ( ti  )exp[-aT&_1  (ti  ) a/2] da 


(4-20) 


£ 
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where  g“[£“Ii(t)&(t)  ]  .  Although  the  integral  of  the  Gaussian 
density  function  does  not  have  a  closed-form  expression,  we 
know  that  the  entire  area  under  a  probability  density 
function  must  have  a  magnitude  of  one, 


U2n)«  |&| ]-i/ «  •  exp ( -fiT R-  1  a/ 2 ) da  -  1 


(4-21) 


Thus,  the  maximum  amplitude  of  the  rate  function,  A(t) ,  can 
be  shown  as: 

A(t)  -  [n/ (tF -to )] • [ ( 2n ) “ I R ( t ) l]-l/a  (4-22) 

where  m  is  the  dimensionality  of  the  vector  a. 

The  mean  noise  arrival  rate  is  defined  as 

fa(t)  -  (t,t)dt  -  (ti  ,r) -k  (4-23) 

where  fa (ti  ,£)  is  not  a  function  of  £  in  this  application 
(fa (t,£)  is  defined  as  being  uniformly  distributed  over  the 
entire  length  of  the  detector) . 

The  data  provided  by  SOFE  is  statistically  reduced  and 
the  results  plotted  by  a  second  program  called  SOFEPL.  The 
statistics  used  to  evaluate  the  performance  of  either  the 
Kalman  or  Meer  filters  are  the  mean  and  standard  deviation 
of  the  error  between  the  truth  and  filter  state  (see 


85 


Equations  (4-2)  and  (4-4)).  Although  the  SOFEPL  program 
provides  an  array  of  statistical  options,  they  are  limited 
to  the  error  statistics  found  between  the  truth  and  filter 
states  (for  filter  performance  evaluation)  and  are 
inadequate  for  evaluating  a  controller  design.  An 
alternative  provided  by  SOFEPL  is  for  the  user  to  specify 
the  desired  error  to  be  evaluated  statistically.  The 
procedure  is  used  to  generate  the  RMS  time  history  as 
described  by  Equation  (4-8) ,  where  the  statistics  are  based 
on  the  error  between  the  true  beam  and  true  target 
positions. 

4.3  The  Beam  and  Target  Parameters 

This  section  reviews  the  beam  filter  and  target  filter 
parameters  since  they  can  play  such  an  important  role  in  the 
design  and  performance  of  the  MMAC.  The  Meer  filter  is 
designed  around  six  parameters.  They  are  as  follows: 

tb  is  the  beam  time  constant. 

g  is  the  square  root  of  the  beam  propagation  noise 
strength,  g  ■  (G2Q)l/a  (Q«l  by  assumption). 

R  is  the  beam  dispersion  measured  as  the  variance  of 
the  Gaussian-shaped  beam  at  the  detector  surface. 

SNR  is  the  signal-to-noise  ratio  (to  be  defined  as  in 
Equation  (4-24) ) . 

D  is  the  MMAE  filter  depth. 

n  is  the  expected  number  of  signal-induced  events 
during  a  simulation  run. 

The  three  parameters  associated  with  the  Gauss-Markov 
acceleration  target  model  are  as  follows: 

tt  is  the  target  time  constant. 


Qt  is  the  target  dynamics  driving  noise  strength. 
Rt  is  the  target  measurement  noise  variance. 

The  signal-to-noise  ratio  is  defined  for  this 
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application  as  the  average  number  of  signal  induced  events 
produced  for  every  noise  event.  The  SNR  can  be  expressed  as 
[18]  : 

SNR  «  [n/(tr-to  )]/[*»  (t,£> -L] 

-  A(2nRM't/[x»  <t,£) -L]  (4-24) 

where  n  is  the  number  of  signal-induced  events  expected 
during  one  simulation  run,  (tf-to)  is  the  duration  of  the 
simulation.  An (t,£)  is  the  noise  arrival  rate  per  length  of 
the  detector  array,  and  L  is  the  length  of  the  detector 
array.  For  this  one-dimensional  application,  the  length  of 
the  detector  is  10  cm. 

The  nominal  values  for  the  beam  and  tracker  filter 
parameters  used  for  the  analyses  are  as  follows: 

The  Beam  Parameters 
T»  *  20  (sec)  D  *  3 

g  *  0.2  (cm2/sec)l/2  R  =  0.5  (cm2) 

SNR  *  20  n  ■  100  (signal  events) 

The  Target  Parameters 

tt  ■  10  (sec)  Qt  *  0.1  (cm*  /sec9)  Rt  *  0.5  (cm2) 
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When  the  the  Meer  filter  is  simplified  and  the  assumptions 
from  Section  2.5  apply,  the  MMAE  filter  depth  is  reduced  to 


one. 

Because  the  fastest  transients  at  the  nominal  condition 
have  a  time  constant  of  ten  seconds  (rather  benign 
dynamics) ,  the  Kalman  filter  and  the  feedback  control  sample 
period  will  be  set  to  one  second.  This  easily  satisfies  the 
Shannon's  sampling  theorem  which  states  that  the  sampling 
rate  should  be  at  least  twice  the  highest  signal  frequency 
content  of  interest.  To  provide  an  extra  margin  of 
insurance,  the  engineering  guide  of  sampling  ten  times 
faster  than  the  highest  frequency  is  used. 

4.4  Performance  Analyses 

Five  different  performance  analyses  are  required  to 
design  and  evaluate  the  MMAC,  which  is  composed  of  a  bank  of 
PI  controllers,  of  which  each  elemental  controller  uses  its 
own  Kalman  filter  to  estimate  the  target  state  and  a  Meer 
filter  to  estimate  the  beam  state.  The  first  two  analyses 
are  the  alternate  controller  cost  weighting  matrices 
analysis  and  the  reduced  Meer  filter  depth  analysis.  The 
results  of  these  analyses  will  allow  us  to  select  the  best 
set  of  cost  weighting  matrices  for  good  on-line  design 
performance  and  adequate  robustness,  and  to  allow  us  to 
justify  using  a  filter  depth  of  one  (i.e. 

D  ■  1) .  The  next  two  analyses  are  the  sensitivity  analysis 
and  the  robustness  analysis.  The  results  from  these 
analyses  are  used  to  develop  the  MMAC.  Last,  another 
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analysis  is  required  to  evaluate  the  operation  of  the  MMAC . 
The  first  four  analyses  are  performed  on  the  PI  tracker 
developed  in  Section  3.3.2.  The  MMAC  analysis  is  performed 
on  the  adaptive  controller  developed  in  Section  3.5.2. 

4.4.1  Evaluating  Alternate  Weighting  Matrices . 

Although  Jamerson  did  some  work  on  selecting  a  suitable  set 
of  controller  cost  weighting  matrices  to  achieve  good 
controller  performance  and  robustness,  an  exhaustive  study 
was  never  accomplished.  Therefore,  a  performance  analysis 
is  conducted  to  determine  the  best  alternate  cost  weighting 
matrices  which  provide  good  performance  at  design  conditions 
and  that  also  enhance  the  MMAC's  robustness.  In  other 
words,  the  goal  is  to  select  the  best  setts)  of  weighting 
matrices  that  best  minimize  the  averaged  RMS  error  (without 
using  excessive  amounts  of  control)  when  the  truth  model 
matches  the  filter  model,  while  simultaneously  providing  at 
least  stability  at  off-design  conditions. 

4.4.2  Evaluating  the  Simplified  Meer  Filter .  The 
second  analysis  evaluates  the  performance  of  the  controller 
and  the  Meer  filter  when  the  filter  depth  is  reduced  from 

D  ■  3  to  D  »  1.  This  will  reduce  the  number  of 
propagating  elemental  Snyder-Fishman  filters  from  eight  to 
two,  and  in  fact  the  filter  can  be  expressed  equivalently 
with  only  a  single  elemental  filter,  but  with  its  gain 
expressed  as  a  function  of  residual  size.  Although  this 
will  slightly  increase  the  average  RMS  error,  it  should 
significantly  decrease  the  computational  loading  (see 
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Section  2.5).  It  makes  on-line  feasibility  much  more 
reasonable . 

4.4.3  Sensitivity  Analysis .  The  purpose  of  the 
sensitivity  analysis  is  to  evaluate  the  performance  of  the 
PI  controller  as  it  is  exposed  to  several  different  "real- 
world"  environments.  Because  each  environment  could  be 
described  by  a  set  of  parameters,  the  sensitivity  analysis 
is  conducted  by  varying  a  single  parameter  in  both  the  truth 
model  and  in  the  controller.  The  results  lend  insight  on 
how  the  controller  reacts  if  it  knows  the  true  environment 
and  the  results  provide  an  absolute  baseline  of  the  best 
possible  performance  that  could  be  expected  from  an  adaptive 
controller.  The  test  is  based  on  the  nominal  parameter 
settings  presented  in  the  previous  sections  and  most  of  the 
parameters  are  evaluated  at  one  order  of  magnitude  above  and 
below  their  nominal  values.  The  test  evaluates  all  the  beam 
and  target  parameters  except  the  filter  depth,  which  is  set 
to  one  (D  *  1) .  The  cost  weighting  matrices  are  based  on 
the  results  in  Chapter  5  and  are  set  to  U=1 ,  Xii=100  and 
Xaa^lO  (refer  to  Equation  (3-62)). 

4.4.4  Robustness  Analysis .  To  complement  the 
sensitivity  analysis,  a  robustness  analysis  is  performed  to 
evaluate  the  controller's  performance  when  an  algorithm- 
assumed  parameter  differs  from  the  "real-world"  environment. 
The  robustness  analysis  is  conducted  under  the  same 
guidelines  as  the  sensitivity  analysis,  but  differs  in  that 
the  parameters  are  changed  in  the  truth  model  without 
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informing  the  controller's  filters.  The  purpose  of  this 
study  and  the  sensitivity  analysis  is  to  identify  those 
parameters  that  are  best  suited  for  on-line  adaptation  by 
the  MMAC,  and  to  identify  any  other  potential  problems  with 
the  PI  controller. 

4.4.5  Evaluating  the  MMAC .  Once  the  MMAC  is 
developed,  it  is  evaluated  against  a  performance  baseline. 
This  baseline  is  defined  as  a  PI  controller  which  receives 
noise-corrupted  measurements,  z(ti),  but  has  access  to  the 
truth  model’s  parameters.  This  baseline  provides  the  best 
possible  performance  the  MMAC  could  achieve  if  it  could 
perfectly  estimate  the  uncertain  parameters.  A  well 
designed  and  tuned  MMAC  should  closely  approximate  this 
baseline.  Also,  the  MMAC  is  compared  to  an  unknowledgeable 
and  non-adaptive  controller  to  see  how  much  better 
performance  the  on-line  parameter  adaptation  yields. 


4 . 5  Summary 

The  purpose  of  this  chapter  has  been  to  explain  the 
tools  and  methods  used  to  develop  and  evaluate  the  MMAC,  and 
to  evaluate  the  performance  of  the  Kalman  and  Meer  filters 
within  the  MMAC  structure.  The  Monte  Carlo  simulation 
provides  the  most  viable  method  of  evaluating  the 
controller's  performance.  SOFE  and  SOFEPL  provide  the  basic 
method  of  computing  the  filter  state  estimation  error 
statistics  and  controller  tracking  error  statistics  which 
allow  us  to  evaluate  the  filters  and  the  controller  designs. 
Four  of  the  performance  analyses,  the  alternate  weighting 
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matrices,  the  reduced  Meer  filter,  the  sensitivity,  and  the 
robustness  analyses,  provide  the  foundation  to  develop  the 
MMAC.  Upon  completion,  the  MMAC  is  developed  and  evaluated 
against  a  baseline  controller,  which  assumes  that  the 
baseline  controller  has  access  to  the  true  parameters  while 
receiving  noise  corrupted  measurements.  The  MMAC  is  also 
compared  to  the  results  of  a  non-adaptive  PI  controller. 

The  results  of  the  five  analyses  are  in  the  next  chapter. 


VS5 


V.  RESULTS  OF  THE  MONTE  CARLO  ANALYSES 


The  results  of  the  five  Monte  Carlo  analyses  discussed 
in  Section  4.4  are  presented  in  this  chapter.  The  first 
analysis  evaluates  seven  different  sets  of  weighting 
matrices  and  found  the  set,  U*l,  Xii*100,  Xos=10,  produces 
the  lowest  steady  state  RMS  tracking  error.  None  of  the 
seven  sets  of  weighting  matrices  tested  cause  the  controller 
to  use  excessive  amounts  of  control.  The  second  analysis 
shows  that  the  controller  can  be  reduced  to  a  filter  depth 
of  one  without  any  measurable  increase  in  RMS  error  between 
the  beam  and  target  positions.  Therefore,  the  rest  of  the 
research  is  conducted  with  the  filter  depth  set  to  one 
(D*l) .  Unfortunately,  an  appreciable  reduction  in  computer 
time  is  not  seen.  This  is  credited  to  the  overall 
inefficiency  of  the  source  code,  and  it  is  assumed  that  a 
program  designed  to  take  advantage  of  a  reduced  filter  depth 
would  require  less  computer  processing  time.  The 
sensitivity  analysis  provides  a  baseline  of  the  best 
possible  controller  performance  provided  the  filter  knows 
the  "real-world”  environment,  while  the  robustness  analysis 
measures  the  ability  of  the  controller  to  cope  with  a 
"real-world"  environment  that  differs  from  the  parameter 
values  assumed  in  the  filter/controller  design.  The 
performance  achieved  when  the  values  of  most  "real  world" 
parameters  are  allowed  to  vary  shows  good  robustness 
characteristics.  The  two  notable  exceptions  are  the  target 


dynamics  driving  noise,  Qr ,  which  the  MMAC  is  designed  to 
handle,  and  the  particle  beam  time  constant,  t» ,  which 
appears  to  be  the  source  of  a  controller  stability  problem. 
Last,  the  MMAC  provides  good  on-line  adaptive  estimation  of 
Qr  .  The  MMAC  significantly  outperforms  the  non-adaptive 
controller,  but  the  controller  responds  slowly  to  a 
decreasing  target  dynamics  noise  strength. 

The  results  of  the  first  four  analyses  are  from 
existing  SOFE  software  used  by  Meer,  Zicker,  Moose  and 
Jamerson  [8,18,20,27],  which  was  modified  for  easier  usage, 
and  corrected  for  three  software  errors.  These  errors  were 
detected  either  by  validation  testing  or  by  reviewing  the 
source  code,  and  the  errors  may  have  affected  the  previous 
results.  The  errors  are  as  follows:  (1)  the  original  SOFE 
source  code,  version  May  1982,  had  a  repeated  line  of  code, 
line  "LJ  »  LJ  +  K"  in  Subroutine  PSQRT.  Subroutine  PSQRT 
computes  the  Cholesky  upper  triangular  square  root  matrix  of 
the  Kalman  filter  covariance  matrix,  £(ti  ) .  This  error 
would  miscalculate  the  upper  triangular  elements  of  the 
Cholesky  square  root  matrix.  (2)  Jamerson's  source  code  had 
an  error  in  the  parameter  list  of  the  call  statement,  "CALL 
RICDSD ( • ) ” .  This  library-obtained  subroutine  calculated  the 
steady-state  solution  to  the  Riccati  equation  used  to 
calculate  the  controller  gains  (see  Equation  3-67) .  This 
error  resulted  in  slightly  higher  controller  gains  than 
would  be  properly  evaluated.  (3)  Zicker' s  source  code, 
which  was  modified  and  used  by  Moose  and  Jamerson, 


calculated  the  true  maximum  amplitude  of  the  rate  function, 
At ,  as  a  function  of  the  filter  beam  dispersion,  Rf  (see 
Equation  (4-22)).  The  true  beam  dispersion,  Rt ,  should  have 
been  used  instead.  This  slightly  affected  the  results  from 
the  previous  robustness  analyses  for  the  beam  dispersion. 

The  results  from  the  MMAC  analysis  are  from  a  corrected 
version  of  the  existing  software  which  is  modified  for 
adaptive  estimation/control.  All  the  subroutines  used  to 
implement  the  MMAC  were  verified  independently  of  SOFE. 


5.1  Evaluating  Alternate  Weighting  Matrices 

The  purpose  of  this  study  is  to  select  the  cost 
weighting  matrix  that  best  minimizes  the  time-averaged  RMS 
error,  RMS« ,  between  the  target  and  the  beam  without  using 
excessive  amounts  of  control.  Therefore,  seven  different 
performance  analyses  have  been  performed  with  a  different 
cost  weighting  matrix  as  defined  by  the  steady-state  cost 
function. 


x*  (ti  ) 

T 

Xi  i  -Xi  i  0  0 

Xa  (ti  ) 

J  -  l  (1/2* 

XT  P  (tl  ) 

-Xii  Xu  0  0 

XT  P  (ti  ) 

q(ti ) 

0  0  Xa  a  0 

q(t» ) 

U(tl  ) 

0  0  0  u 

■ 

u(tl  ) 

The  results  as  depicted  by  Table  5-1  and  5-2  show  that  the 
best  cost  weighting  matrix  is  composed  of  U*l,  Xi  i -100  and 
Xa  a "10 .  A  closer  analysis  shows  that  a  10:1  ratio  of  Xi i  to 
Xa  a  is  the  important  relationship  that  future  designs  should 


TABLE  5-1 


The  RMS  Errors  from  the  Cost  Weighting  Matrices  Analysis 


u 

Xi  1 

Xa  s 

Average 

RMS* 

Minimum 
RMS*  (ti  ) 

Maximum 
RMS.  (tt  ) 

% 

1 

1 

1 

2.536 

2.258 

2.881 

+3.6% 

1 

10 

10 

2.639 

2.314 

2.918 

+7.8% 

1 

10 

100 

2.929 

2.526 

3.251 

+19.7% 

1 

100 

10 

2.448 

2.165 

2.737 

- 

1 

100 

100 

2.676 

2.336 

2.963 

+9.3% 

1 

100 

1000 

2.956 

2.546 

3.286 

+20.7% 

1 

1000 

100 

2.450 

2.166 

2.740 

+0.1% 

TABLE  5-2 

The  Applied  Control  During  the  Second  Monte  Carlo  Run 


U 

X!  1 

Xa  a 

RMS  of 
u(tt  ) 

Minimum 

U(tl  ) 

Maximum 
u(ti  ) 

* 

1 

1 

1 

18.363 

-50.896 

18.046 

3 

1 

10 

10 

18.413 

-51.153 

18.688 

4 

1 

10 

100 

18.586 

-52.456 

20.231 

6 

1 

100 

10 

18.330 

-51.276 

17.451 

1 

1 

100 

100 

18.440 

-51.365 

18.753 

5 

1 

100 

1000 

18.608 

-52.478 

20.470 

7 

1 

1000 

100 

18.333 

-51.283 

17.740 

2 

Note:  *  is  the  order  of  performance  with  #1  being  the 

best . 
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consider.  This  cost  weighting  matrix  allows  the  controller 
to  operate  with  the  lowest  average  RMS  error  among  the  seven 
sets  evaluated.  Although  none  of  the  cost  weighting 
matrices  allow  the  controller  to  use  significantly  large 
amounts  of  control  in  relation  to  the  severity  of  the  target 
dynamics,  the  best  cost  weighting  matrix  provides  the  least 
oscillatory  application  of  control  as  measured  by  the 
applied  control  RMS  statistic  (see  Table  5-2) .  These 
findings  are  consistent  with  the  theory  presented  in  Section 
3-4.  As  we  increase  the  pseudo-integrator  cost  weighting 
element,  Xao,  with  respect  to  the  position  states  cost 
weighting  element,  Xi i ,  the  controller  will  begin  to 
overcompensate  large  differences  in  the  target  and  beam 
positions.  This  will  cause  the  controller  to  be  more 
oscillatory  in  the  application  of  control  and  it  will  lead 
to  larger  RMS  errors.  Also,  as  the  cost  weighting  of 
control,  U,  approaches  the  values  of  the  other  cost 
elements,  the  control  will  become  more  expensive  to  apply, 
and  this  will  tend  to  raise  the  average  RMS  error.  A  plot 
of  the  RMS  error  statistic  is  plotted  in  Figure  5-1.  Except 
for  analyzing  the  pole  placement  of  the  controller  (refer  to 
Chapter  6) ,  nothing  was  done  to  enhance  the  robustness  of 
the  controller  in  any  explicit  manner. 

The  relationship  between  the  amount  of  applied  control 
and  the  relative  error  between  the  target  and  the  controlled 
beam  is  shown  in  Figure  5-2.  This  plot  is  from  the  second 
simulation  run  out  of  the  Monte  Carlo  analysis  and  shows  how 


the  applied  control  reacts  to  the  target.  In  this  example, 
the  target  is  accelerating  across  the  detector  array  of  the 
tracking  particle  beam,  which  causes  the  magnitude  of 
applied  control  to  increase  continuously.  The  plot  starts 
at  ti “50  seconds  because  the  first  50  seconds  is  being  used 
to  let  the  system's  transient  response  attenuate.  The 
discontinuities  in  the  relative  error  are  do  to  the  Kalman 
filter  updates. 

The  previous  concerns  [8]  over  the  effect  of  a  large 
initial  overshoot  caused  by  the  beam  position  being  offset 
from  the  target  are  important  because  the  initial  transient 
response  can  saturate  the  actuators  and  be  very 
destablizing.  They  were  ignored  only  because  the  emphasis 
is  on  the  steady  state  operation  of  the  estimator/ 
controller,  and  it  is  suggested  that  during  target 
acquisition  a  different  cost  weighting  matrix  or  filter 
tuning  parameter  be  used  which  is  more  forgiving  of  large 
initialization  errors. 

5.2  Evaluating  Simplified  Meer  Filter 

The  purpose  of  this  analysis  is  to  evaluate  the 
performance  of  the  controller  when  the  depth  of  the  Meer 
filter  is  reduced  from  D»3  to  D*l.  The  goal  is  to  reduce 
the  computational  loading  of  the  Meer  filter  significantly 
at  the  cost  of  slightly  degrading  the  performance  of  the 
controller.  The  insight  was  provided  by  Zicker  [27]  who 
found  that  when  the  Meer  filter  was  reduced  from  a  depth  of 
3  to  1,  the  average  RMS  error  increased  by  only  0.004  cm. 
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The  results  from  the  controller’s  performance  analyses  show 
that  the  controller  is  even  more  tolerant  to  the  changes  in 
filter  depth.  As  the  data  in  Table  5-3  suggest,  there  is  no 
detectable  degradation  of  the  controller's  performance  when 
the  depth  of  the  Meer  filter  is  reduced  to  one. 
Unfortunately,  the  expected  savings  in  CPU  time  did  not 
appear . 


TABLE  5-3 

Results  of  the  Reduced  Meer  Filter  Analysis 


Average 

Minimum 

Maximum 

CPU 

Depth 

U 

Xi  1 

Xa  s 

RMS. 

RMS.  (ti  ) 

RMS.  (ti  ) 

Time 

3 

1 

100 

10 

2.44795 

2.16530 

2.736H 

511.5 

1 

2.44795 

2.16531 

2.73671 

521.6 

3 

1 

100 

100 

2.67559 

2.3357£ 

2.96271 

500.5 

1 

2.67559 

2.33575 

2.96271 

497.0 

Note:  CPU  Time  is  for  the  entire  simulation,  which  includes 

data  reduction  and  plotting;  units  are  in  seconds. 
Underlining  of  digits  is  used  to  accentuate  the 
differences  in  computed  statistics;  units  are  in  cm. 


Because  these  findings  did  raise  some  suspicion,  the 
analysis  was  supported  by  an  in-depth  review  of  the  source 
code.  The  review  indicated  that  the  code  was  correct,  and 
that  the  reason  there  was  no  apparent  reduction  in  CPU  time 
was  due  to  the  general  inefficiency  of  the  source  code  which 
did  not  take  advantage  of  the  reduced  filter  depth. 
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One  of  the  reasons  that  contributes  to  the  outstanding 
performance  of  the  controller  with  the  reduced  Meer  filter 
is  the  simple  beam  dynamics  used  by  both  the  plant  and 
filter.  A  more  realistic  plant  might  have  had  a  much  more 
complicated  transfer  function  which  could  lead  to  greater, 
more  discernible  errors  when  the  Meer  filter's  depth  is 
reduced.  In  other  words,  a  depth  of  one  might  be  sufficient 
for  a  first-order  Gauss-Markov  scalar  position  process,  but 
a  greater  depth  might  be  required  if  a  higher-order  Gauss- 
Markov  model  were  a  significantly  better  model.  Another 
reason  is  that  the  high  SNR  (SNR*20)  provides  few  noise- 
induced  events  to  test  the  Meer  filter.  Combining  the 
effects  of  a  high  SNR  with  the  effects  of  a  tight  beam 
dispersion  (R=0.5  cm2)  and  a  long  detector  array  (L-10  cm), 
results  in  the  fact  that  most  of  the  noise  events  should  be 
easily  detected  as  being  induced  by  noise  rather  than 
signal,  and  the  Meer  filter  does  not  require  the  extra 
filter  depth  to  accomplish  this  adaptive  decision-making. 

5.3  Sensitivity  Analysis 

The  purpose  of  the  sensitivity  analysis  is  to  evaluate 
a  "fully  parameter-knowledgeable"  controller,  in  which  both 
the  Kalman  and  Meer  filters  (as  well  as  the  controller  gain 
computations)  have  access  to  the  true  parameters.  The 
results  provide  the  best  possible  performance  that  can  be 
expected  of  the  controller  under  various  "real  world" 
conditions,  and  the  results  define  the  baseline  for  the 
robustness  analysis.  Throughout  the  sensitivity  and 


robustness  (Section  5.4)  analyses,  the  Meer  and  Kalman 
filter  parameters  are  set  to  the  the  nominal  conditions 
defined  by  Table  5-4,  and  then  each  nreal  world"  parameter 
is  evaluated  individually  through  variations  plus  or  minus 
one  order  of  magnitude.  In  these  sensitivity  studies,  the 
corresponding  filter/controller-assumed  parameter  is  varied 
accordingly,  while  in  the  later  robustness  studies,  the 
filter/controller-assumed  parameter  is  left  unchanged.  The 
sensitivity  results  from  the  beam  parameters  are  presented 
first. 

TABLE  5-4 

Nominal  Conditions  for  the  Sensitivity  and  Robustness 

Analyses 


Beam 

Parameters  (see 

Section  4.3) 

T»  ■ 

20  (sec) 

D  *  1 

0  - 

0.2  (cm2 /sec) 1 ' 2 

R  »  0. 5 

(cm2  ) 

SNR  « 

20 

n  ■  100 

(signal  events) 

Target  Parameters  (see  Section  4.3) 

Tt  * 

10  (sec) 

Rt  *  0.5 

(cm2  ) 

Or  - 

0.1  (cm2 /sec) 

Non-Zero  Elements  in  the  Cost  Weighting  Matrix 
(see  Section  3.4  and  Equation  3-62) 

U  -  1  Xu  -  100  Xoo  *  10 


5.3.1  Beam  Time 


-  ia_.  The  beam  time  constant 


represents  the  relative  speed  with  which  the  particle  beam 
dynamics  can  react  to  the  applied  control.  The  effect  tb 
has  on  the  plant  dynamics  can  be  seen  by  means  of : 


B  B 

Gx  (s)  -  -  -  -  (5-2) 

s  -  pi  s  +  (1/t«  ) 

where  Gi (s)  is  the  plant  transfer  function.  The  effects  of 
increasing  or  decreasing  tb  has  the  effect  of  moving  the 
pole  closer  to  or  away  from  the  origin  of  the  s-plane.  The 
effects  of  moving  the  pole  can  be  evaluated  by  looking  at 
the  plant's  settling  time  for  a  step  input, 

ts  -  3 . 912/ I  pi  I  *  3.912*Tb  (5-3) 

where  ts  is  defined  as  the  time  it  takes  the  exponentially 
time-correlated  plant  dynamics  to  reach  98  percent  of  the 
final  steady  state  value  [3] .  As  tb  is  decreased,  the  pole 
moves  farther  from  the  origin  of  the  s-plane,  the  beam 
dynamics  become  quicker,  and  it  should  be  easier  for  the 
particle  beam  to  track  the  target.  This  should  result  in  a 
smaller  average  RMS  error  and  this  is  confirmed  by  the 
results  in  Table  5-5.  But,  contrary  to  what  one  might 
expect,  the  large  changes  in  tb  do  not  result  in  large 
changes  in  the  tracker  error.  This  insensitivity  is  due  to 
the  fact  that  the  system's  frequency  response  (as  can  be 
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observed  in  a  Bode  plot)  is  not  only  a  function  of  the 


s 


% 


9 


a 

% 


i  i 


I  l 


isKsswr 


corner  frequency,  wc  (uc»1/ti  )  ,  but  also  it  is  a  function  of 
the  loq  magnitude  of  the  open  loop  gain.  Km  (Km*Btb ) ,  and 
the  effects  of  changing  tb  are  cancelled  by  changes  in  Km . 
The  filter's  tb  represents  the  time  constant  in  the  filter's 
model  of  the  plant,  whereas  the  true  tb  is  the  actual  time 
constant  of  the  particle  beam  plant.  In  a  sensitivity 
analysis,  the  true  and  filter  parameter  values  are  equal. 


TABLE  5-5 

The  Sensitivity  Analysis  Results  for  tb 


tb  (seconds) 

Average 

Filter 

True 

RMS. 

% 

Comments 

2 

2 

2.431 

i 

o 

• 

<#> 

Quicker  beam  dynamics 

20 

20 

2.448 

- 

Nominal  conditions 

200 

200 

2.456 

+0.3% 

Slower  beam  dynamics 

5.3.2  Beam  Propagation  Noise  Strength  -  g»-  The  beam 
propagation  noise  strength  indicates  the  confidence  we  have 
in  the  target  model  being  correct  and  the  magnitude  of  the 
random  fluctuations  that  the  beam  position  actually 
undergoes.  As  we  increase  the  strength  of  the  noise,  g* ,  we 
are  increasing  the  magnitude  of  the  random  beam  position 
fluctuations,  and  the  Meer  filter  must  depend  more  on  the 
measurement  update.  This  can  be  demonstrated  by  analyzing 
this  relationship  with  an  elemental  Snyder-Fishman  filter 
covariance  propagation  equation,  and  gain  equation: 


r 


P(tl*l- ) 


,tt  )P(ti*  )  +  Qd  (ti  ) 


(5-4) 


K(ti)  -  P(ti- )/(P(ti- )  +  R)  (5-5) 

where  the  closed-form  solution  of  Qd  (ti )  is 

Qd (ti )  =  g2 tb [1-exp (-2At/Ti )] /2  (5-6) 

where  At  is  the  varying  sample  time  between  the  arrival  of 
events.  As  g2  increases,  so  does  P(ti»i~)  increase; 
therefore,  K(ti )  will  increase  toward  a  value  of  one,  and 
the  filter  will  begin  to  rely  more  heavily  on  the 
measurement  update. 

With  a  sensitivity  analysis,  we  assume  that  the  filter 
g2  matches  a  hypothetical  true  g* ,  which  is  used  in  the  true 
state  dynamics  model, 


§ 

i 

£ 


fe 


XB  t  ( ti  ♦  1  )  »  Ob  t  ( ti  ♦  1  ,  tt  )  Xb  t  ( ti  )  +  Bd  ( ti  )  U  ( ti  )  +  Wd  t  ( ti  ) 

(5-7) 

where  Wdt(tt)  is  a  zero-mean,  white  Gaussian,  discrete-time , 
stochastic  process  with  a  variance  as  given  in  Equation 
(5-6) .  As  we  increase  the  true  g2  (and  the  filter  g2 
correspondingly) ,  the  magnitude  of  the  random  fluctuations 
beam  undergoes  increases,  and  the  filter  must  depend  more 
heavily  on  the  measurements  which  come  from  the  stochastic 
truth  model.  Therefore,  as  the  true  g2  increases  and  the 
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actual  RMS  value  of  excursions  of  the  beam  increase,  so 


should  the  average  RMS  error.  The  results  in  Table  5-6 
confirm  this. 


TABLE  5-6 

The  Sensitivity  Analysis  Results  for  g2 


g  (cm2 /sec)1 / 2 
Filter  True 

Average 

RMS. 

% 

Comments 

0.02 

0.02 

2.419 

-1.2% 

K(ti )  decreases 

0.2 

0.2 

2.448 

- 

Nominal  conditions 

2.0 

2.0 

3.743 

+53.9% 

K(ti )  increases 

5.3.3  Beam  Dispersion  -  R  The  particle  beam  is  said 
to  have  a  Gaussian-shaped  dispersion  as  measured  with 
respect  to  the  surface  of  the  photodetector.  Because  the 
beam  is  assumed  to  have  Gaussian  distribution  (see  Equation 
2-1) ,  it  can  be  represented  statistically  by  the  true  mean, 
t  (t) ,  and  the  true  "variance",  Rt  (t) .  As  Rt  (t)  increases, 
the  signal-induced  events  spread  over  a  larger  area  of  the 
detector,  making  it  more  difficult  to  estimate  the  center  of 
the  beam,  and  forcing  the  filter  to  rely  more  heavily  on  its 
internal  dynamics  model.  Although  Table  5-7  supports  this 
concept,  the  filter  is  relatively  insensitive  to  equal 
changes  in  both  the  true  and  filter  R's. 


TABLE  5-7 
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The  Sensitivity  Analysis  Results  for  R 


R  (cm2 ) 

Filter  True 

Average 

RMS* 

% 

Comments 

0.05 

0.05 

2.437 

<#> 

• 

o 

i 

Depend  on  z> (ti ) 

0.5 

0.5 

2.448 

- 

Nominal  conditions 

5.0 

5.0 

2.475 

+1.1% 

Depend  on  filter 

5.3.4  Sional-to-Noise  Ratio  -  SNR  The  signal-to- 
noise  ratio  is  the  ratio  of  the  number  the  signal  events  to 
the  number  of  noise  events  that  arrive  over  a  period  of 
time.  Within  the  context  of  this  simulation,  the  filter  SNR 
is  used  to  estimate  the  noise  rate  parameter  (the  SNR  is 
defined  in  Equation  (4-24), 

XNf(t,t)  *  [n/  ( tr  -to  )  ]  /  [SNRf  •  L]  (5-8) 

The  true  SNR  is  used  to  generate  the  the  actual  noise 
parameter  rate  used  to  generate  the  noise  events.  Thus, 
when  the  SNR  is  decreased,  we  expect  to  see  the  increase  of 
noise  events  to  drive  up  the  time-averaged  RMS  error,  as 
depicted  in  Table  5-8.  The  table  shows  a  discrepancy  in 
that  the  time-averaged  RMS  error  shows  a  slight  increase 


A 


when  the  SNR  is  increased.  After  reviewing  the  results  from 
the  robustness  analysis  (see  Section  5.4.4),  one  can  see 
that  the  controller  is  highly  insensitive  to  changes  in  SNR, 
and  that  the  small  deviations  of  time-averaged  RMS  error  is 
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within  the  sample  statistical  error.  This  is  also  supported 
by  the  fact  that  the  different  noise  arrival  rates  will 
cause  different  random-number  generated  signal-  and  noise- 
induced  event  patterns  (i.e.,  the  same  random  number 
generator  seed  was  used  throughout  the  study,  and  changes  in 
signal  or  noise  parameter  rates  will  alter  the  simulation 
pattern  of  signal-  and  noise-induced  events  by  requiring 
different  numbers  of  calls  to  the  random  number  generator) . 


TABLE  5-8 

The  Sensitivity  Analysis  Results  for  SNR 


SNR 

Filter 

True 

Average 

RMS. 

% 

Comments 

2 

2 

2.469 

+0.9% 

1/3  noise  events 

20 

20 

2.448 

- 

Nominal  conditions 

200 

200 

2.466 

+0.7% 

Nearly  noise  free 

5.3.5  Expected  Number  o£  Signal-Induced  Events  -  a 
The  expected  number  of  signal  events  is  used  to  calculate 
the  filter's  maximum  amplitude  of  the  rate  function,  as  by 
Equation  (4-22) , 

Af  -  nt /[<2flRf  )‘/i  •  (tr-to  )]  (5-9) 


which,  in  turn  is  used  by  the  filter  to  calculate  the 
conditional  probability  that  a  signal-induced  event  did 
occur  (see  Equations  (2-1)  and  (2-27) .  The  true  n  is  used 
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to  calculate  the  signal  arrival  rate,  xs  (t)  .  As  both  the 
true  and  filter  n  are  decreased  by  equal  amounts,  the  filter 
must  rely  more  heavily  on  the  filter  dynamics  model  because 
the  mean  Meer  filter  sample  period  is  increased.  Therefore, 
we  should  expect  the  time-average  RMS  error  to  increase  as 
shown  in  Table  5-9.  The  RMS  error  statistics  may  be 
slightly  affected  by  differences  in  the  random-number 
generated  signal-  and  noise-induced  event  patterns  caused  by 
a  change  in  the  signal  parameter  rate.  The  n=1000 
performance  run  was  not  accomplished  because  the  expected 
results  did  not  warrant  the  high  computer  costs.  More 
signal  events  should  slightly  increase  the  filter's 
accuracy.  The  n=10  performance  run  should  have  produced  an 
unstable  system  (according  to  findings  in  Chapter  6)  because 
the  mean  sample  period  of  10  seconds  would  place  the  poles 
of  the  discrete-time  controller  far  outside  the  unit  circle, 
but  the  perfectly  matched  filter  time  constants  prevented  us 
from  seeing  the  instability. 


TABLE  5-9 

The  Sensitivity  Analysis  Results  for  n 


5.3.6  Target  Time  Constant  -  Tr_  The  target 
acceleration  time  constant  is  used  in  the  target  shaping 
filter  which  simulates  the  target  dynamics.  The  shaping 
filter  is  a  linear,  time-invariant  system  driven  by  a 
stationary  white  Gaussian  noise  process  that  has  a  mean  of 
zero  and  a  variance  of  Qt  .  The  target  dynamics  position 
process  is  modeled  as  the  double  integral  of  exponentially 
time-correlated  acceleration;  a  third-order  Markov  process 
that  is  the  output  of  a  third  order  shaping  filter  which  is 
shown  in  Figure  5-3.  By  definition,  the  white  noise  process 


1/tt 


Figure  5-3.  Target  Shaping  Filter 


has  a  power  spectral  density  (PSD)  magnitude  of  Qt  and  an 
infinitely  wide  bandwidth  as  shown  in 

Figure  5-4a.  When  this  white  noise  is  passed  through  a 
first-order  lag  filter,  as  is  done  in  the  target  shaping 


Figure  5-4.  PSD  Functions:  a)  White  Noise,  b)  Target 

Acceleration 


filter,  it  becomes  the  acceleration  of  the  target,  and  has  a 
power  spectral  density  function  of 

2ot  2  /tt  Qt 

PSDa(w)  =  -  =  -  (5-10) 

w2  +  (1/tt  ) 2  w2  +  (  1/tt  ) 2 

where  the  variance  (mean  squared  value)  of  the  target 
acceleration  is  defined  as 


Ot2  *  (RMSa)2  *  HQtTt 


(5-11) 


A  plot  of  the  acceleration  power  spectral  density  function 
is  shown  in  Figure  5-4b.  The  effects  of  reducing  tt  would 
increase  the  bandwidth  of  the  target  acceleration  and  lead 
to  faster  dynamics,  but  with  a  much  smaller  amplitude  as 
shown  in  Figure  5-5.  The  overall  results  of  a  small  tt  is  a 
lower  RMSa  (assuming  Qt  is  held  constant) .  Therefore, 
smaller  tt  should  lead  to  smaller  time-averaged  RMS  errors 
as  shown  in  Table  5-10. 

TABLE  5-10 

The  Sensitivity  Analysis  Results  for  Tt 


tt 

Filter 

(sec) 

True 

Average 

RMS. 

% 

Comments 

1.0 

1.0 

1.325 

-45.9% 

ot  2 "0 . 05 ,  PSD=0.1 

2.5 

2.5 

1.882 

-23.1% 

Or  2 =0.125,  PSD=0 . 625 

5.0 

5.0 

2.225 

-9.1% 

Ox  2  *0 . 25 ,  PSD=2.5 

10.0 

10.0 

2.448 

- 

Nominal:  ot2=)4,  PSD=10 

100.0 

100.0 

2.685 

+9.7% 

Or  2  *5 . 0 ,  PSD  =1000 

PSD  is  calculated  for  PSD(w=0) ,  i.e.,  the  low 
frequency  asymptotic  value  at  w»0. 


5.3.7  Target  Dynamics  Noise  Strength  -  Qt  The  target 
dynamics  noise  strength  is  the  strength  of  the  white 
Gaussian  noise  used  to  drive  the  shaping  filter.  Although 
it  does  not  effect  the  bandwidth  of  the  target  acceleration, 
it  directly  affects  the  amplitude  of  the  power  spectral 
density  function  as  can  be  seen  by  setting  w  to  zero  in 
Equation  (5-10) : 

PSD ( w=»0)  *  QtTt*  (5-12) 

Therefore,  an  increase  in  Qt  will  increase  the  RMS*  (see 
Equation  (5-11))  and  should  result  in  an  increase  in  RMS 
error.  Table  5-11  confirms  this.  The  true  Qt  is  the  actual 
white  noise  strength  used  to  drive  the  shaping  filter,  where 
the  filter  Qt  is  used  by  the  Kalman  filter  in  the  target 
state  estimate  and  covariance  propagation  equations. 


TABLE  5-11 

The  Sensitivity  Analysis  Results  for  Qt 


Qt  (cm*  /sec8  ) 
Filter  True 

Average 

RMSe 

% 

Comments 

0.01 

0.01 

1.461 

-40.3% 

ot  * “0 . 05 ,  PSD-1 

0.1 

0.1 

2.448 

- 

Nominal:  ot*«H,  PSD-10 

1.0 

1.0 

4.644 

+89.7% 

Ot  8  =5. 0 ,  PSD-100 

Note:  PSD  is  calculated  for  PSD(w«0) . 


5.3.8  Evaluating  tt  and  Or  with  the  Target 
Acceleration  Power  Spectral  Density  Maximum  Set  to  10 . 0 
Although  it  appears  that  the  filter  is  much  more  sensitive 
to  Qt ,  one  more  analysis  is  required  to  confirm  this.  For 
this  analysis,  the  amplitude  of  the  low  frequency  asymptote 
of  PSD  (see  Figure  5-5)  is  held  constant,  and  tt  is  varied. 
The  results  in  Table  5-12  show  that,  for  targets  with  the 
same  low  frequency  acceleration  PSD  characteristics,  those 
which  can  exhibit  more  violent  maneuvering  (i.e.,  the  tt  is 
smaller,  or  the  acceleration  process  bandwidth  1/tt  is 
larger)  are  much  more  difficult  to  track. 


TABLE  5-12 

The  Sensitivity  Analysis  Results  for  PSD 


PSD  * 

10.0  (cm2 /Hz)  at 

u*0 

Filter 

True 

Average 

tt 

Qt 

tt 

Qt 

RMS. 

% 

1.0 

10.0 

1.0 

10.0 

5.384 

+119.9% 

10.0 

0.1 

10.0 

0.1 

2.448 

- 

100.0 

0.001 

100.0 

0.001 

1.089 

-55.5% 

Note:  PSD  is  calculated  for  PSD(w*0). 

5.3.9  Target  Measurement  Noise  Variance  -  Rj_  The 
target  measurement  noise  variance  indicates  our  confidence 
in  the  measurement  model.  As  the  measurements  becomes  more 
severely  corrupted,  the  measurements  assume  a  wider  and 


flatter  Gaussian  distribution  centered  around  the  true 


target  position,  and  the  variance,  which  is  represented  by 
the  true  Rt ,  increases.  The  filter  Rt  is  the  filter's  value 
representing  the  actual  Rt .  As  the  filter  Rt  increases,  the 
Kalman  filter  gain  decreases  and  the  filter  must  rely  more 
heavily  on  its  internal  dynamics  model,  and  tracking  is  more 
difficult  with  less  precise  sensors.  Therefore,  an  increase 
in  Rt  should  result  in  a  higher  time-averaged  RMS  error. 
Table  5-13  confirms  this. 


TABLE  5-13 

The  Sensitivity  Analysis  Results  for  Rt 


Rt  (cm2  ) 
Filter 

True 

Average 

RMS. 

% 

Comments 

0.05 

0.05 

1.518 

-38.0% 

Relies  on  measurements 
which  are  precise 

0.5 

0.5 

2.448 

- 

Nominal  conditions 

5.0 

5.0 

4.473 

+82.7% 

Relies  on  filter’s 
model  and  measurements 
are  imprecise 

5.4  Robustness  Analysis 

The  purpose  of  the  robustness  analysis  is  to  evaluate 
the  controller's  performance  when  the  embedded  Kalman  or 
Meer  filters  mismodels  the  "real  world".  Initially,  all  the 
parameters  are  set  to  the  nominal  conditions,  and  then  the 


true  parameters  are  evaluated  one  at  a  time  at  values  above 
and  below  the  nominal  condition  to  simulate  the  effects  of 
the  filters  mismodeling  the  "real  world".  The  results  from 


the  robustness  analysis  are  evaluated  against  a  baseline 
defined  by  the  results  of  the  sensitivity  analysis. 

Tracking  performance  is  highly  insensitive  to  the  changes  in 
most  of  the  true  beam  parameters,  and  by  definition,  this  is 
a  demonstration  of  excellent  robustness  characteristics  (at 
the  chosen  nominal  conditions) .  The  notable  exception  is 
the  beam  time  constant,  for  which  the  results  indicate  a 
possible  stability  problem.  Some  poor  robustness 
characteristics  were  demonstrated  for  all  three  target 
parameters,  of  which,  the  target  dynamics  noise  strength  was 
selected  as  the  parameter  best  suited  for  on-line  adaptive 
estimation. 

5.4.1  Beam  Time  Constant  -  Tj_  Originally,  the  true  tb 
was  to  be  evaluated  at  plus  and  minus  one  order  of  magnitude 
of  the  nominal  filter  tb  so  that  the  results  would 
correspond  to  the  sensitivity  analysis,  but  severe 
computational  difficulties  led  to  run-time  computer  failures 
before  the  simulations  could  be  completed.  To  be  able  to 
collect  useful  data,  the  true  tb  was  evaluated  at  plus  five 
percent  and  minus  four  percent  of  the  nominal  filter  tb .  In 
all,  five  robustness  analyses  were  conducted  with  true  tb 
equal  to  19.2,  19.6,  20.0  (nominal  conditions),  20.4,  and 
21.0  seconds.  The  results  of  the  analyses  are  presented  in 
Table  5-14.  The  respective  plots  of  RMS  errors  between  the 
target  and  the  beam  position,  and  the  mean  and  standard 
deviation  of  the  error  between  the  true  and  filter  beam 
state  are  shown  in  Figures  15-6  through  15-15. 


The  results  indicate  a  severe  stability  robustness 
problem.  When  the  true  tb  is  simulated  at  the  values  19.2 
or  21.0  seconds,  the  results  in  the  both  table  and  in  the 
plots  indicate  that  the  controller's  performance  is 
suffering  from  a  strong  instability.  Growing  RMS  values  and 
standard  deviations  that  increase  with  the  severity  of  the 
mismodeling  strongly  support  the  concept  that  a  stability 
problem  exists.  These  findings  led  to  a  sensitivity 
evaluation  of  the  controller  gains  and  a  stability  analysis 
of  the  PI  controller. 


TABLE  5-14 


The 

Robustness 

Analysis 

Results  for  tb 

tb 

Filter 

(sec) 

True 

Average 

RMS. 

Minimum 
RMS.  (ti  ) 

Maximum 
RMS.  (ti  ) 

Tuned 

RMS. 

% 

20 

19.2 

80.080* 

4.555 

338.690* 

«2 . 448 

*  * 

20 

19.6 

4.779 

2.920 

7.030 

s2.448 

+195% 

20 

20.0 

2.448 

2.165 

2.737 

2.448 

- 

20 

20.4 

4.687 

2.949 

7.062 

*2.448 

+191% 

20 

21.0 

96.719* 

5.301 

364.960* 

=2 .448 

*  * 

Note:  %  This  is  the  percentage  of  difference  between  the 
average  RMS  error  and  the  tuned  RMS  error,  which 
is  based  on  or  approximated  from  the  results  of 
the  sensitivity  analysis  (i.e.,  tb  t  =ti  f ) . 

*  These  values  and  the  results  from  the  Meer  filter 
covariance  analysis  (see  Figure  5-7)  indicate  that 
the  controller  frequently  became  unstable  before 
the  completion  of  the  run  (before  100  seconds) . 


Extremely  large  values. 


i 


Figure  5-7.  Mean  and  ±  Standard  Deviation  Envelope  of 
Error  between  the  True  and  Filter  Position  States: 
ti f  *20 . 0  (sec),  Ti t  *19 . 2  (sec) 
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Figure  5-11.  Mean  and  ±  Standard  Deviation  Envelope  of  the 
Error  between  the  True  and  Filter  Position  States: 
T*f-20.0  (sec),  Tit«20.0  (sec) 
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Figure  5-12.  RMS  Error  Between  the  Target  and  Bean 
Position:  Tif-20.0  (sec),  Tat-20. 4  (sec) 
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Figure  5-13.  Mean  and  ±  Standard  Deviation  Envelope  of  the 
Error  between  the  True  and  Filter  Position  States: 
t*  f  *20 . 0  (sec),  Ti t “20 . 4  (sec) 
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Figure  5-15.  Mean  and  ±  Standard  Deviation  Envelope  of  the 
Error  between  the  True  and  Filter  Position  States: 
Tif-20.0  (sec),  Tit-21.0  (sec) 
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The  reason  for  the  controller  gam  sensitivity 
evaluation  is  to  see  if  the  gains  are  sensitive  to  changes 
in  tb .  If  the  controller  gains  are  found  to  be  highly 
sensitive  to  changes  in  tb ,  it  is  possible  that  the  poles  of 
the  full-state  feedback  system  are  being  driven  outside  the 
unit  circle,  because  of  the  gains,  Gc* ,  rather  than  due  to 
the  effects  of  the  Meer  filter  being  in  the  loop.  This 
requires  that  the  controller  gain  equations.  Equations 
(3-65)  through  (3-70)  be  solved  for  different  values  of  tb . 
The  results  in  Table  5-15  suggest  that  the  controller  gains 
are  relativity  insensitive  to  changes  in  tb . 

The  results  of  the  stability  analysis  are  in  Chapter  6, 
indicating  that  the  stability  robustness  problem  is  induced 
by  mismodeling  of  the  real  world  in  the  Meer  filter  dynamics 
model.  Since  this  is  a  totally  different  form  of  analysis, 
it  is  separated  into  a  different  chapter.  Because  we  do  not 
know  the  cause  of  the  instability  (at  this  point  in  time) , 
the  beam  time  constant  would  not  be  suitable  for  on-line 
adaptive  estimation. 

TABLE  5-15 

Sensitivity  of  the  Controller  Gains  to  Changes  in  tb 


Tb 

Gen*  Gc  12*  Gc  l  3*  Gc  l  4  * 

Gca* 

E 

19.0 

1.2413  -1.2918  -1.0261  -0.48H 

0.2747 

1.0760 

20.0 

1.2421  -1.2921  -1.0252  -0.4880 

0.2744 

1.0746 

21.0 

1.2431  -1.2906  -1.0240  -0.4874 

0.2740 

1.0733 

controller  is  rather  insensitive  to  mismodeling  of  the  beam 
propagation  noise  strength,  as  is  shown  in  Table  5-16. 
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TABLE  5-16 

The  Robustness  Analysis  Results  for  g2 


g2  (cm2 /sec) 
Filter  True 


0.02 

0.2 

2.0 


Average 

Tuned 

RMS. 

RMS. 

% 

2.427 

2.419 

+0.3% 

2.448 

2.448 

- 

4.032 

3.743 

+7.7% 

5.4.3  Beam  Dispersion  -  R  The  controller  is  rather 
insensitive  to  mismodeling  of  the  beam  dispersion  parameter 
as  can  be  seen  from  the  results  in  Table  5-17. 


TABLE  5-17 

The  Robustness  Analysis  Results  for  R 


(cm2  ) 

ter  True 

Average 

RMS. 

Tuned 

RMS. 

% 

.5  0.05 

2.439 

2.437 

+0.1% 

.5  0.5 

2.448 

2.448 

- 

.5  5.0 

2.549 

2.475 

+3.0% 

5.4.4  Signal  to  Noise  Ratio  -  SNR  The  results  in 


Table  5-18  indicate  that  the  controller  is  completely 
insensitive  to  SNR  parameter  mismodeling. 


TABLE  5-18 

The  Robustness  Analysis  Results  for  SNR 


SNR 

Filter 

True 

Average 

RMS. 

Tuned 

RMS. 

% 

20 

2 

2.496 

2.496 

0.0% 

20 

20 

2.448 

2.448 

- 

20 

200 

2.466 

2.466 

0.0% 

5.4.5  Expected  Number  of  Signal-Induced  Events  -  n 
The  controller  is  slightly  sensitive  to  mismodeling  of  the 
number  of  signal  events,  as  the  results  in  Table  5-19 
indicate.  Note  concerning  the  data  —  the  tuned  (baseline) 
RMS  error  for  nt  *200  is  only  an  approximation  and  no  data 
was  collected  at  that  value  during  the  sensitivity  analysis 
(see  Table  5-9) .  The  effect  is  depicted  in  Table  5-19  as 
approximations  in  the  "Tuned  RMS. "  column.  Also,  changes  in 
the  signal  arrival  rate  will  cause  different  random-number 
generated  signal-  and  noise-induced  event  patterns,  which 
could  have  contributed  to  the  statistical  differences  in  the 
time-averaged  RMS  errors. 


TABLE  5-19 


The  Robustness  Analysis  Results  for  n 


n  (signal 
Filter 

events) 

True 

Average 

RMS. 

Tuned 

RMS. 

% 

100 

10 

2.522 

s2.540 

<#> 

r~~ 

• 

o 

i 

100 

100 

2.448 

2.448 

- 

100 

200 

2.468 

s2.448 

+0.8% 

5.4.6  Taroet  Time  Constant  -  u_  The  controller 
appears  to  be  highly  sensitive  to  target  time  constant 
mismodeling  when  true  target  time  constant  is  less  than  the 
filter  time  constant  (see  Table  5-20) .  This  is  as 
anticipated,  since  a  small  true  tt  is  a  harsher  target  to 
track  than  one  with  a  long  true  tt .  This  lack  of  robustness 
indicates  that  the  controller  underreacts  to  the  much 
quicker  target  even  though  the  target  RMS  acceleration  has 
decreased  substantially  (refer  to  Table  5-10  or  Equation 
(5-11)  to  see  the  relationship  between  tt  and  RMS 
acceleration).  Meanwhile,  a  much  slower  target  dynamics 
characteristic  in  the  real  world  allows  the  controller  to 
track  the  target  easily.  The  improvement  over  the  tuned 
(baseline)  value  is  not  easily  explained.  It  is  doubtful 
the  entire  1.3  percent  difference  can  be  explained  as  a 
statistical  error  since  both  the  baseline  sensitivity 
analysis  and  the  robustness  analysis  used  the  same  signal- 
and  noise-induced  event  pattern. 


5.4.8  Bv 


aluatina  Tt  and  Or  with  the  Target 
Acceleration  Power  Spectral  Density  Maximum  Set  to  10.0 
In  this  robustness  study,  the  true  parameters  tt  and  Qr  are 
varied  while  the  maximum  PSD  was  held  to  a  constant  value  of 
10.0  cm  as  defined  by  Equation  (5-12).  The  results  in  Table 
5-23  confirmed  the  overall  controller  sensitivity  to  Qt ,  but 
the  findings  are  not  as  extreme  as  the  results  from  the 
previous  two  sections  might  have  suggested. 


TABLE  5-22 

The  Robustness  Analysis  Results  for  Constant  PSD 


PSD  »  10. 
Filter 

Tt  Qt 

0  ( cm* /Hz ) 
True 

Tt  Qt 

Average 

RMS. 

Tuned 

RMS. 

% 

10 

1 

10.0 

6.810 

5.384 

+26.5% 

10 

1 

10 

0.1 

2.448 

2.448 

- 

10 

100 

0.001 

1.995 

1.089 

+83.2% 

Note:  PSD  is  calculated  for  PSDCw^O) . 


5.4.9  Target  Measurement  Noise  Variance  -  Rl  Table 
5-22  shows  that  the  controller  is  most  sensitive  to  Rt 
mismodeling  when  the  measurement  of  the  target  state  is  much 
more  severely  corrupted  by  noise  than  the  filter  assumes. 
This  is  expected.  The  importance  of  this  analysis  lies  with 
the  target's  ability  to  deceive  or  corrupt  tracking  by  way 
of  jamming  the  controller's  sensors. 


TABLE  5-23 


The  Robustness  Analysis  Results  for  Rt 


Rt 

Filter 

(cm2  ) 

True 

Average 

RMS* 

Tuned 

RMS* 

% 

0.1 

0.01 

1.623 

1.518 

+6.9% 

0.1 

0.1 

2.448 

2.448 

- 

0.1 

1.0 

6.324 

4.473 

+41.4% 

5.5  Evaluating  the  Multiple  Model  Adaptive  Controller 

This  section  presents  the  results  obtained  during  the 
evaluation  of  the  MMAC  design  described  in  Section  3.5. 

The  adaptive  controller  is  composed  of  three  elemental  PI 
controllers,  each  tuned  to  a  prespecified  uncertain 
parameter  value,  Qt ,  the  target  model  white  noise  strength. 
The  first  and  third  elemental  controllers  are  based  on  the 
target's  upper  and  lower  limits  of  Qt ,  where  Qt i  represents 
the  white  noise  strength  of  a  benign  trajectory  (Qt 1=0.01 
cm2 /sec8 ) ,  and  Qt3  is  the  maximum  white  noise  strength 
associated  with  evasive  maneuvering  (Qts=1.0  cm2 /sec8 ) .  Thi 
middle  elemental  controller  is  based  on  an  intermediate  Qt 
equal  to  0.07  cm2 /sec8 ,  as  opposed  to  the  previous  nominal 
value  of  0.1  cm2 /sec® .  This  value  of  dynamics  target  noise 
strength  was  selected  from  a  series  of  trial  tests,  each 
composed  of  20  simulation  runs.  It  was  found  that  as  Qt2 
approached  either  the  upper  or  lower  limits,  the  adaptive 
controller  had  difficulty  distinguishing  Qt2  from  the  other 
Qt's,  and  the  controller  error  increased. 


A  successful  adaptive  controller  must  meet  two 
important  criteria.  First,  the  adaptive  controller  must 
improve  robustness  significantly  to  warrant  the  additional 
computational  load  of  an  adaptive  structure.  This  can  be 
evaluated  by  comparing  the  Monte  Carlo  analysis  results  of 
the  MMAC  against  a  baseline  performance,  the  sensitivity 
analysis  results,  and  the  robustness  analysis  results.  The 
baseline  performance  is  the  best  possible  performance  that 
can  be  expected  of  the  adaptive  controller.  The  baseline 
performance  is  done  by  informing  the  adaptive  controller  of 
the  true  Qt  and  setting  the  probabilistic  weighting  of  the 
"correct"  elemental  controller  to  one.  The  baseline 
performance  in  this  evaluation  used  the  filter-estimated 
states,  because  the  filter  states  were  used  to  generate  the 
results  found  in  Sections  5.1  through  5.4.  (An  alternative 
method  could  assume  an  "ultimate”  baseline  controller  to 
have  access  to  the  truth  states,  but  this  is  an  artificial 
assumption  that  yields  an  unrealistic  baseline  of 
performance . ) 

The  Monte  Carlo  analysis  for  the  actual  adaptive 
controller  and  baseline  controller  uses  the  following  test 
format:  from  t=0.0  to  t=50.0  seconds,  the  true  Qt  is  set  to 

0.01,  and  the  first  40.0  seconds  are  used  to  let  the  system 
transients  die  out.  At  t=50.0  seconds,  the  true  Qt  is 
changed  to  1.00  and  the  target  enters  its  evasive 
maneuvering  phase  of  the  test.  At  t»80.0  seconds,  the  true 
Qt  is  reset  to  0.01,  and  the  target  resumes  a  benign 


trajectory.  Data  is  taken  over  three  periods  of  time.  The 
overall  RMS  error  is  calculated  by  averaging  from  t=40.0  to 
100.0  seconds,  while  the  controller's  ability  to  adapt  to 
the  benign  trajectory  is  measured  by  averaging  from  t=40.0 
to  50.0  seconds,  and  the  controller's  ability  to  adapt  to 
the  evasive  maneuvering  is  measured  by  averaging  from  t=70.0 
to  80.0  seconds.  The  speed  and  effectiveness  of  the 
adaptation  is  indicated  from  the  the  time  history  of  the  RMS 
values,  which  can  be  observed  in  Figures  5-16  and  5-17. 

The  results  in  Table  5-24  show  that  the  adaptive 
controller  does  an  excellent  job  of  estimating  the  true  Qt 
representing  the  evasive  maneuvering,  where  performance  is 
most  critical.  Large  errors  here  could  result  in  the  loss 
of  lock  on  the  target  (i.e.,  the  target  moves  off  the 
detector  array) .  Also,  the  adaptive  controller  does  a  good 
job  of  adapting  to  the  Qt  appi  opriate  for  a  return  to  a 
benign  trajectory,  although  the  performance  is  not  as  good. 
The  data  provided  by  the  sensitivity  analysis  provides  a 
baseline  of  the  best  performance  the  PI  controller  could 
attain  if  it  had  access  to  the  true  Qt .  The  differences 
between  this  baseline  and  the  MMAC  baseline  demonstrate  the 
degradation  in  tracker  performance  due  to  using  on-line 
adaptation  and  selecting  only  three  discretized  values  of  Qt 
for  the  MMAC  structure.  The  data  provided  by  the  robustness 
analysis  shows  what  the  results  would  have  been  if  we  did 
not  use  on-line  adaptive  estimation  of  Qt . 
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5-16.  MMAC 
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the  Beam 
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Qt  i*0.01 
Qt  2  =0 . 07 
Qt  3  *1 . 00 
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Average 

RMS. 

MMAC 
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RMS. 
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tivity: 

Average 

RMS. 

Robustness : 

Average  RMS. 

Qt  =0 . 01  Qt  =0 . 1  Qt=1.0 

Overall 
40< t<100 

3.139 

3.546 

- 

- 

Qt  *0.01 
40< t<50 

1.634 

1.736 

1.461 

2.034  3.340 

Qt  =1 . 00 
70<t<80 


4.657 


4.643 


4.644 


7.300  4.998 


smallest  magnitude  of  variance  residual,  An  (ti ) .  If  this 
happens,  the  controller  probabilistic  weighting  will  be 
erroneous,  because  the  magnitudes  of  An (ti  )  are  independent 
of  both  the  residuals  and  the  correctness  of  the  different 
elemental  controllers  [1,133. 

Although  the  actual  residual  values  of  the  filters 
within  the  adaptive  controller  structure  were  not  collected, 
the  probabilistic  weightings  were  recorded.  The  plots  in 
Figure  5-18  are  from  the  first  two  Monte  Carlo  runs  and 
suggest  the  adaptive  controller  had  little  difficulty  in 
distinguishing  the  "correct''  elemental  filter/controller 
from  the  "mismodeled”  f ilter/controllers . 

The  second  important  criterion  that  should  be  used  when 
evaluating  an  adaptive  tracker  is  how  quickly  the  tracker 
can  detect  and  respond  to  changes  in  target  dynamics  noise 
strength.  By  definition,  parameters  are  more  slowly  varying 
(if  they  vary  at  all)  than  the  states  [13] .  Unfortunately, 
the  true  Qt  is  controlled  by  the  target,  and  although  it 
varies  slowly  when  it  is  compared  to  the  target  states,  the 
target  can  quickly  change  Qt  with  respect  to  the  one  second 
controller  sample  period  used  in  this  study.  (This  is 
unrealistically  long,  but  has  been  chosen  to  be  constant 
with  the  previous  feasibility  studies  [8,20,27].) 

Therefore,  it  is  important  to  see  how  fast  the  MMAC  can 
adapt  to  changes  in  Qt .  The  data  is  taken  from  a  single 
simulation  run  that  lasted  2,000  seconds,  with  the  initial 
Qt  set  to  0.01  cm* /sec8 .  At  65  seconds,  and  at  every  100 


seconds  thereafter,  the  true  Qt  is  increased  to  1.0 
cm2 /sec9 ,  and  at  90  seconds  and  at  every  100  seconds 
thereafter,  the  true  Qt  resumes  0.01  cm2 /sec2 . 

The  results  in  Table  5-25  indicate  that  the  adaptive 
controller  takes  approximately  five  seconds  to  detect  a 
change  in  Qr ,  and  it  is  somewhat  slow  in  adapting  to  the  new 
Qt .  In  addition,  it  takes  the  adaptive  controller  much 
longer  to  react  to  a  decrease  in  Qt  than  it  takes  to  react 
to  an  increase  in  Qt .  When  you  go  from  benign  to  a  harsh 
trajectory,  [rit2/Ak]  of  the  benign  filter  k  becomes  very 
large  while  the  corresponding  term  in  the  filter  based  on  a 
harsh  trajectory  model  is  small,  so  the  probabilities 
calculated  in  Equation  (3-79}  causes  the  probability 
weighting  to  shift  to  the  latter  filter.  However,  when  the 
actual  trajectory  switches  back  to  a  benign  trajectory,  the 
value  of  [rk2/Ak]  in  the  mismatched  filter  is  not 
extraordinarily  large  compared  to  that  of  the  matched 
filter,  so  the  probability  weighting  shift  is  not  as  rapid 
[11,15].  This  effect  can  be  observed  in  Figures  5-16  and  5- 
17.  A  possible  solution  to  quicken  the  adaptive 
controller's  reaction  to  changes  in  Qt  is  to  increase  both 
the  controller  and  Kalman  filter  sample  rates.  The  other 
problem  of  the  filter/controller  slowly  adapting  to  a 
decreasing  Qt  is  inherent  to  the  design,  and  has  been 
previously  noted  in  other  studies  of  multiple  model  adaptive 
algorithms  [11,15]. 
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TABLE  5-25 


Adaptive  Filter  Time  Response  to  Increasing 
and  Decreasing  Qt 


A.  Qt  increases  from  0.01  to 

»-* 

• 

o 

o 

• 

Sample 

Time  responses  in  seconds 

Statistics 

to 

t» 

ts 

tT 

£  (mean) 

5.30 

8.45 

6.90 

13.75 

o  (std  dev) 

2.27 

8.11 

6.70 

7.04 

B.  Qt  decreases  from  1.00  to 

0.01. 

Sample 

Time  response  in  seconds 

Statistics 

to 

t«  * 

ts  * 

tT  * 

£  (mean) 

5.05 

31.5 

20.5 

36.5 

o  (std  dev) 

2.06 

15.3 

25.8 

15.2 

Note:  The  sample  population  is  n*20,  except  when  denoted 

by  *,  then  n*19  (i.e.,  the  2000  second  run  from  which 
this  data  is  tabulated,  terminated  before  the 
information  on  t« ,  t s  and  tT  was  available) .  The 
sampled  data  is  in  Appendix  D. 

Definitions: 

to  -  the  delay  time  is  the  time  it  takes  the  adaptive 
controller  to  react  first  to  a  new  Qt . 
ts  -  the  settling  time  is  the  time  required  by  the 

"correct"  elemental  filter  to  reach  a  pic  >0.5, 
and  stay  above  it  thereafter.  It  is  measured 
from  the  time  the  adaptive  controller  first 
reacts  to  the  new  Qt . 

tt  -  the  rise  time  is  the  time  it  takes  the  "correct" 
elemental  filter  to  surpass  a  pic  of  0.75. 

Note:  This  is  within  10  percent  of  the  upper 
bound.  It  is  measured  from  the  time  the 
adaptive  controller  first  reacts  to  the  new  Qt . 
tr  -  the  transition  time  is  the  time  it  takes  the 

"correct"  elemental  filter  to  detect  and  react 
to  a  change  in  Qt  .  tT  =  to  +  t* 
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Five  different  Monte  Carlo  analyses  were  accomplished. 
The  first  four  evaluated  the  non-adaptive  PI  controller 
while  the  fifth  examined  the  MMAC .  It  was  found  that  the  PI 
controller  performed  best  with  a  10:1  ratio  between  the  cost 
weighting  elements  Xi 1  and  X99,  and  that  the 
estimator/controller  suffered  no  detrimental  effects  when 
the  Meer  filter's  depth  was  reduced  from  three  to  one.  The 
performance  was  found  to  be  the  most  sensitive  to  equal 
changes  in  Qt  (the  target  dynamics  driving  noise  strength) 
and  Rt  (the  measurement  corruption  noise  variance)  in  both 
the  truth  model  and  the  filters,  while  tb  (the  beam  time 
constant)  and  Qt  were  the  parameters  that  yielded  the  most 
significant  robustness  difficulties.  The  severe  lack  of 
robustness  characteristics  with  respect  to  tb  appears  to  be 
a  stability  problem  which  will  be  discussed  in  Chapter  6. 

The  target  driving  noise  strength  was  selected  for  on-line 
adaptive  estimation,  and  the  adaptive  controller  performed 


VI.  PI  CONTROLLER  STABILITY  ANALYSIS 


The  purpose  of  this  chapter  is  to  evaluate  how  the  Meer 
and  Snyder-Fishman  filters  affect  the  stability  of  the  PI 
controller.  The  motivation  behind  this  analysis  is  to 
isolate  the  stability  robustness  problem  mentioned  in 
Section  5.4.1,  and  if  possible,  identify  the  reasons  the  PI 
controller  lacks  robustness  with  respect  to  the  beam  time 
constant.  The  source  of  the  instability  is  believed  to  lie 
within  one  of  three  areas:  with  the  mismodeling  of  the 
filter's  state  transition  matrix  (O)  ,  within  the  inherent 
structure  of  the  Meer  filter,  or  with  the  varying  sample 
period  which  is  the  product  of  the  Poisson  process  for 
measurement  event  times. 

The  chapter  will  start  by  evaluating  the  stability  of 
the  PI  full-state  feedback  controller  structure  that  uses 
the  steady-state  gains  found  in  Table  3-1.  The  generic  PI 
controller  structure  is  time-invariant  and  linear; 
therefore,  showing  that  the  poles  of  the  system  lie  inside 
the  unit  circle  (  -domain)  will  prove  that  the  closed  loop 
system  with  a  full-state  feedback  controller  structure  is 
asymptotically  stable  in  the  global  sense  (a  form  of  zero- 
input  stability) .  Then,  the  Meer  filter  is  cascaded  with 
the  controller  transfer  function,  and  the  stability  is 
evaluated  as  functions  of  the  Meer  filter  sample  rate  and 
the  filter  gain.  Because  the  sample  rate  of  the  Meer  filter 
is  not  fixed,  the  system  is  no  longer  time-invariant,  such 


a  simple  stability  evaluation  is  not  possible.  A  rigorous 
proof  of  global  stability  for  stochastic  control  systems 
that  use  either  the  Meer  or  Snyder-Fishman  filters  in 
conjunction  with  Kalman  filters  does  not  exist.  Instead,  an 
ad  hoc  method  is  employed  where  the  effects  of  different 
sample  rates  are  evaluated  as  being  time-invariant,  and  the 
zero-input  stability  analysis  is  used  to  evaluate  the  pole 
movement  within  the  stochastic  controller.  Next,  a 
stability  demonstration  is  used  to  verify  that  an  analogous 
Kalman  filter  for  this  application  is  stable  and  inferences 
are  drawn  to  the  possible  sources  of  instability  within 
Snyder-Fishman  filter.  Last,  we  need  to  analyze  the  effects 
the  Meer  filter  structure  has  on  the  Snyder-Fishman  filter 
equations,  and  thus  the  effects  the  Meer  filter  structure 
has  on  the  PI  controller.  The  emphasis  in  this  chapter  is 
not  to  provide  the  rigorous  conditions  for  global  stability, 
but  to  analyze  the  possible  causes  of  the  controller's 
instability.  In  other  words,  we  are  looking  for  answers  to 
the  findings  in  Section  5.4.1. 


6.1  Derivation  ol 


Controller  Transfer  Functions 


In  a  linear,  time-invariant  system,  we  find  that  the 
stability  is  determined  by  the  location  of  the  poles  in  the 
system  transfer  function.  A  system  is  said  to  be 
asymptotically  stable  if  and  only  if  all  the  poles  lie 


within  the  stable  region.  For  the  continuous  system,  the 
stable  region  is  defined  as  the  left-hand  plane  of  the 
Laplace  s-domain  plot.  Analogously  for  the  discrete  system 


I 


model,  the  stable  region  is  defined  as  the  area  inside  of 
the  unit  circle  of  a  z-domain  plot.  In  other  words,  a 
discrete,  linear,  time-invariant  system  model,  such  as  the 
closed  loop  system  with  deterministic  full-state  feedback  PI 
controller,  is  asymptotically  stable  if  and  only  if  all  the 
poles  of  the  discrete  closed  loop  transfer  function  lie 
within  the  unit  circle  (i.e.,  have  a  magnitude  less  than 
one)  . 

The  first  step  is  to  analyze  the  stability  of  the 
deterministic  PI  controller  loop  shown  in  Figure  6-1.  The 
transfer  function  is  formed  by  solving  the  pseudo-integrator 
equation. 


q(z)  *  -  [yc  (z)  -  yr(z)] 

z  -  1 


feedback  control  equation, 


u  ( z )  =  -Gc  i  yc  ( z )  -  Gc  z  q  ( z )  +  yr  ( z ) 


(6-1) 


(6-2) 


and  the  plant  equation  in  terms  of  yc(z)/yr(z).  The  plant 
equation  is  found  by  taking  the  impulse  transform  of  a  zero- 
order  hold  attached  to  the  plant  dynamics  model, 


Yc  (s)  »  Gz  o ■  ( s )  P  { s )  U*  (s) 


(6-3) 


where  Gz  ot (s)  is  the  transfer  function  (in  the  Laplace 
domain)  representing  a  zero-order  hold  and  P(s)  is  the  plant 
transfer  function 


Gz  o  ■ ( s )  =  [1  -  exp(ats)]/s 


(6-4) 


1 

P (s)  =  - 

s  -  l/r»t 


and  then  transforming  it  into  the  z-domain, 


(6-5) 


tb  t  [l-exp(-at/TB t )] 

yc  (z)  =  GzotP(z)u(z)  - - u(z)  (6-6) 

z  -  exp ( —at/ tb  t ) 


where 

GzoiP(z)  -  Zl[l  -  exp (-st) /s] [1/ (s  +  l/TBt)]l  (6-7) 

The  deterministic  PI  controller  transfer  function  is 

Ye  (Z)  Tb  t  [1-exp (-At/TB t ) ] [z (Gc  z  +1) -1] 

yr (z)  [z-exp(-at/TB t ) ] [z-lj 

+  Tb  t  [1-exp (-At/TB »  ) ]  [z (Gc 1 +Gc  a ) -Gc 1 ] 

(6-8) 

Assuming  we  want  to  know  the  location  of  the  poles  for  the 
nominal  conditions  (at-1.0  seconds,  TBt*=20.0  seconds, 
Gci*~1.2422,  and  Gc a  * *0 . 27436) ,  Equation  (6-8)  becomes 
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y«  (z)  . 7654 [z  -  .7847] 

-  -  -  (6-9) 

yr  (z)  z2  -  .  4719z  +  .2605 

The  poles  of  the  closed  loop  system  with  deterministic  PI 
full-state  feedback  controller  transfer  function  are  located 
at  pi *0. 7982  and  p* =-0.3263  which  are  well  inside  the  unit 
circle,  and  therefore  the  deterministic  controller  is 
asymptotically  stable. 

The  next  step  is  to  derive  the  stochastic  PI  controller 
loop  transfer  function  with  Meer  filter  in  the  loop.  For 
simplicity,  the  Kalman  filter,  which  is  used  to  estimate  the 
reference  variable,  will  be  removed  from  the  overall 
controller  structure,  and  its  stability  will  be  verified  in 
Section  6.3.  This  can  be  done  because  the  Kalman  filter 
transfer  function  is  not  included  in  the  feedback  loop  of 
the  regulator  associated  with  the  tracker.  Therefore, 
instead  of  analyzing  the  stability  of  the  PI  tracker,  we 
will  analyze  the  stochastic  PI  regulator,  and  this 
simplification  will  not  effect  the  pole/zero  movement  of 
interest  [14] . 

The  stochastic  controller  loop  transfer  function  is 
derived  in  the  same  manner  as  was  the  deterministic 
controller  transfer  function,  except  that  an  approximation 
to  a  Meer/Snyder-Fishman  filter  is  included  in  the  feedback 
loop  (see  Figure  6-2) .  Note  that  the  approximation  being 
made  is  that  a  fixed  sample  period  is  being  used  rather  than 
taking  measurements  at  times  dictated  by  a  Poisson  process. 
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Figure  6-2.  Stochastic  PI  Controller 


Thus,  the  approximation  results  in  a  Kalman  filter  structure 
as  shown  in  Figure  6-2.  Such  an  approximation  is  made  only 
in  an  attempt  to  discern  basic  stability  robustness 
characteristics  of  the  loop,  possibly  as  a  function  of  a 
chosen  fixed  sample  period,  to  gain  insights  into  the  actual 
behavior  of  a  Meer  or  Snyder-Fishman  filter  in  the  loop. 
Therefore,  the  transfer  function  is  formed  by  solving  the 
equations  for  the  approximated  Meer (or  Snyder-Fishman) 
filter, 

zK ( z ) yc  ( z )  +  Bdu(z)  tl-K(z)  ] 

£♦  (z)  =  -  (6-10) 

z  -  Ob [1-K ( z ) ] 

the  pseudo-integrator, 

z 

q ( z )  =  -  [£+  (z)  -  yr(z)]  (6-11) 

z  -  1 

the  feedback  control, 

u(z)  *  -Gci£Mz)  -  Gc  2  q  ( z )  +  yr  (z)  (6-12) 

and  the  plant  model.  Equation  (6-6),  in  terms  of 
yc(z)/yr(z).  The  term  Bd  is  the  discrete-time  control  input 
matrix  and  can  be  calculated  from  Equation  (3-60) .  The 
stochastic  PI  controller  closed  loop  transfer  function  is: 


yc  (z)  [r*t  ( 1-exp (- At / tb t ) ) ] [z (Gc 2 +1) -1] [z-$b  (l-K(z) ) ] 


yr  (z)  z3  +  z2  [A+D-C]  +  Z  [B-AC-E]  -  BC 

(6-13) 

where 

A*  [l-K(z)]  [Bd  (Gci+Gc2  )-0B]-l 
B  =  [l-K(z)  ]  [Ob  -Bd  Gc  i  ] 

C  =  exp (-At/TB  t ) 

D  *  Tb t K (z) [Gc i +Gc 2 ] [1-exp (-At/TB t ) ] 

E  »  Tb t K(z) [Gc l ] [1-exp (-At/TB t ) ] 


6.2  The  Zero-Input  Stability  Analysis 

The  purpose  of  the  zero-input  stability  analysis  is  to 
see  what  conditions  will  drive  the  stochastic  PI  controller 
unstable.  Two  variables  of  the  Meer  filter  will  be 
analyzed.  First,  the  stochastic  filter  will  be  set  to  a 
steady-state  condition  with  a  "fixed"  sample  period  of  one 
second.  This  is  to  simulate  a  rather  well  behaved  Poisson 
distributed  sample  process  in  which  the  time  "between" 
events  has  closely  approximated  the  mean  signal  arrival  rate 
of  one  second  between  events.  The  reasoning  is  that  we  need 
to  define  some  sort  of  baseline  to  be  able  the  calculate  a 
filter  gain.  Then,  we  want  to  vary  the  time  to  the  "next" 
signal-induced  event  (for  the  varying  sample  period 
analysis,  we  will  discount  the  effects  that  noise-induced 
events  have  on  the  simplified  Meer  filter  gain)  and  observe 


the  e££ects  that  different  sample  periods  have  on  the 
controller's  pole  locations.  In  other  words,  we  want  to  see 
in  an  approximate  manner  how  the  Poisson  distributed  sample 
rate  will  affect  the  controller's  stability.  Theoretically, 
as  the  time  to  the  "next”  event  is  increased,  one  or  more 
poles  will  migrate  beyond  the  limits  of  the  unit  circle,  and 
the  controller  is  forced  to  operate  in  the  unstable  region 
for  the  duration  of  that  sample  period  [28] . 

Second,  the  controller's  stability  will  be  evaluated 
with  respect  to  changes  in  the  Meer  filter  gain,  K(ti ) . 
Unlike  Kalman  or  Snyder-Fishman  filters  gains,  the 
simplified  Meer  filter  gain  (assuming  a  filter  depth  of  one) 
is  a  function  of  the  filter's  residuals  (see  Section  2.5); 
therefore,  the  Meer  filter  is  susceptible  to  large  swings  in 
filter  gain  due  to  noise-induced  events.  As  it  will  be 
seen,  this  can  have  a  significant  effect  on  the  stability  of 
a  controller. 

The  effects  that  a  varying  sample  rate  and  changing 
filter  gains  have  on  controller  stability  can  be  evaluated 
by  solving  for  the  stochastic  PI  controller  closed  loop 
transfer  function,  Equation  (6-13) ,  for  a  range  of  different 
sample  periods  and  filter  gains.  This  is  done  by  solving 
the  state  transition  matrix,  <Pb  ,  for  the  "next"  sample 
period.  The  control  input  matrix,  Bd ,  is  calculated  from 
Equation  (3-60)  for  the  controller  sample  rate  of  one  Hertz 
(i.e.,  a  sample  period  of  at*l  second).  For  the  evaluation 
of  the  varying  sample  rate,  the  filter  gain  is  calculated 
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from  the  Snyder-Fishman  covariance  equations  (see  Equations 
(2-11)  and  (2-14),  where  P(ti-i*)  is  calculated  for  the 
steady-state  condition  defined  by  a  "fixed"  sample  period  of 
one  second  and  P(ti-)  is  calculated  for  the  duration  of  the 
"next"  sample  period.  The  filter  gain  is  calculated  from 
P(ti-)  using  the  Snyder-Fishman  filter  gain  equation. 
Equation  (2-15) .  Changes  in  the  filter  gain  are  evaluated 
by  solving  the  stochastic  PI  controller  transfer  function 
for  a  range  of  different  filter  gains  at  different  sample 
periods . 

6.2.1  The  Effects  of  a  Varying  the  Sample  Period  It 
was  found  that  the  PI  feedback  controller  (using  Meer  filter 
to  estimate  the  beam  states,  and  operating  under  nominal 
conditions)  is  driven  unstable  for  sample  periods  that  are 
longer  than  3.368  seconds.  This  can  be  observed  from  the 
stochastic  controller  root  locus  plotted  in  Figure  6-3. 
Although  sample  periods  longer  than  3.368  seconds  occurs 
only  3.45  percent  of  the  time  (as  calculated  from  Equation 
(4-17),  and  assuming  is *1.0  signal-induced  event  per 
second) ,  it  means  that  the  filter  is  operating  in  the 
unstable  region  for  15.05  percent  of  the  time  (as  calculated 
by  integrating  over  the  Poisson  density  function  from 
ti*3.368  to  ti  *!■•). 

The  effect  of  mismodeling  t»  was  also  studied,  and  the 
results  showed  (see  Appendix  C)  that  the  controller  becomes 
unstable  at  shorter  sample  periods  when  the  true  t«  is 
decreased  from  the  Meer  filter-assumed  value. 


Figure  £-3.  Root  Locus  of  ths  Stochastic  PI  Controller; 
Varying  the  Heer  Filter  Saaple  Rate 


6.2.2  The  Effects  of  a  Varying  Stochastic  Filter  Gain 
The  analysis  was  done  with  the  sample  period  set  to  1.0  and 
3.5  seconds,  and  the  gain  was  varied  from  0.0  to  1.0.  That 
is,  the  filter  gain  in  this  analysis  was  calculated  across 
its  entire  range,  and  the  filter  gain  was  not  based  on  a 
specific  filter  beam  dynamic  driving  noise  strength,  g2 ,  nor 
a  measurement  noise  variance,  R.  The  results  for  a  sample 
period  of  one  second  demonstrated  that  the  controller  is 
stable  for  the  entire  range  of  the  filter  gain  (see  Table 
6-1).  This  is  expected,  for  we  effectively  evaluated  the 
controller  as  if  it  were  using  a  Kalman  filter  instead  of  a 
Meer  filter.  The  controller  is  third  order,  but  the  filter 
employs  a  pole-zero  cancellation  (during  ideally  modeled 
conditions  and  not  true  for  robustness) ,  and  the  effects  of 
the  filter  within  the  controller  structure  is  never  felt. 
This  is  not  the  case  when  the  controller's  sample  period  is 
3.5  seconds  long.  The  controller  is  only  stable  for  values 
of  K(ti )  that  are  less  than  0.27.  This  is  a  relatively 
small  filter  gain,  which  indicates  that  the  filter  assumes 
that  its  internal  dynamics  model  is  relatively  accurate 
during  that  sample  period.  The  results  suggest  that  the 
stability  of  the  controller  is  highly  sensitive  to  the 
magnitude  of  the  filter  dynamic  driving  noise  strength,  g2 , 
which  is  used  to  calculate  the  filter  gain.  The  higher  the 
g2 ,  the  more  susceptible  the  stochastic  controller  is  to 
being  driven  unstable  by  long  sample  periods  (particularly 


if  the  filter-assumed  model  does  not  match  the  real  world 


exactly,  which  is  always  the  case) . 


TABLE  6-1 

Stochastic  Filter  Pole/Zero  Placement  as  a  Function  of  K(ti ) 


Pole  locations  (nominal  condition) 


Filter 

Gain 

at 

Pi 

■  1.0  seconds 

P2  P3 

at 

Pi 

=  3.5  seconds 

P2  P3 

wwm 

.7982 

-.3263 

.0000 

.8183 

-3.8484 

.0000 

H 

.7982 

-.3263 

.0951 

.8183 

-3.4365 

.0111 

|  1 

.7982 

-.3263 

.2378 

.8183 

-2.8248 

.0338 

.7982 

-.3263 

.4756 

.8183 

-1.8375 

.1032 

0.25 

.7982 

-.3263 

.7134 

.8183 

-.9711 

.2949 

0.10 

.7982 

-.3263 

.8561 

.8183 

-.6075 

.5657 

0.00 

.7982 

-.3263 

.9512 

.8183 

-.4560 

.8395 

Zero  locations 

(Note : 

the  Z2  /p3 

cancellation  at 

at=1.0) 

Filter 


1 . 0  seconds 
Z2 


3.5  seconds 

Z2 


.0000 

.0951 

.2378 

.4756 

.7134 

.8561 

.9512 


7847 

7847 

7847 

7847 

,7847 


.0000 

.0839 

.2099 

.4197 

.6296 

.7555 

.8395 


6.3  The  Stability  of  the  Kalman  and  Snvder-Fishman  Filter 
Equations 

This  section  applies  the  second  Lyapunov  stability 
theory  to  the  homogeneous  portion  of  the  Kalman  filter  to 
show  that  it  meets  the  zero-input  stability  criteria.  Then, 
the  same  stability  criteria  will  be  applied  approximately  to 


m 
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the  Snyder-Fishman  filter  and  although  the  stochastic 
particle  beam  model  is  technically  stable,  an  observability 
problem  is  apparent.  The  stability  criteria  to  be  used  in 
this  section  presents  a  set  of  non-restrictive  conditions 
for  which  the  Kalman  filter  algorithm  is  asymptotically 
stable  in  the  global  sense.  Asymptotic  global  stability  is 
the  strongest  form  of  stability,  which  means  that,  given  any 
initial  condition  within  a  defined  region,  its  state  plane 
trajectory  will  always  converge  to  an  equilibrium  point  [3] . 
If  that  system  is  linear  and  asymptotically  stable,  there 
will  be  only  one  equilibrium  point.  This  is  an  important 
concept,  because  for  linear  systems,  the  sufficient 
conditions  of  asymptotic  stability  are  indistinguishable 
from  conditions  for  bounded  input-bounded  output  (BIBO) 
stability  [12] .  When  BIBO  stability  criteria  are  applied  to 
the  filter,  it  means  that,  for  a  bounded  input  into  the 
filter,  the  output  (state  estimate)  will  be  bounded.  If  the 
region  of  the  initial  conditions  that  meet  the  stability 
criteria  is  expanded  to  include  the  entire  state  space,  then 
that  system  is  said  to  be  stable  in  the  global  sense.  If 
the  system  model  upon  which  the  Kalman  filter  is  based  is 
stochastically  controllable  with  respect  to  the  points  of 
entry  of  the  dynamics  driving  noise  w*  ( • ,  • )  ,  and 
stochastically  observable  from  the  points  of  exit  of  the 
measurements  from  the  system  model,  then  the  filter  is 
asymptotically  stable  in  the  global  sense  [12,14]. 


Therefore,  a  Kalman  filter  as  defined  by  the  target 


k  • 
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model  in  Section  3.1,  is  said  to  be  stochastically 
controllable  if  there  exist  positive  numbers  a  and  3, 
0<a<(S<~,  and  a  positive  integer  N,  such  that,  for  all 
i*N, 


al  s  z  e(ti  ,  tj  )£<i  (tj-i  )fid  (tj-i  )GdT  (tj-i  )eT(ti  ,tj  )  s  pi 

J  «  t  -  N  ♦  t 

(6-14) 

and  is  said  to  be  stochastically  observable  if  there  exist 
positive  numbers  a  and  3,  0<a<p<-,  and  a  positive  integer  N, 
such  that,  for  all  iiN 


alS  I  OMtj  ,ti  )&Mtj  )&-»  (tj  )H(tj  )$(tj  ,ti  j  i  PI  (6-15) 


Because  Equations  (6-14)  and  (6-15)  are  bounded  both  above 
and  below,  this  is  a  stronger  condition  than  just  using  the 
positive  definiteness  of  the  controllability  and 
observability  gramians.  If  the  target  model  upon  which  the 
Kalman  filter  is  based  is  stochastically  controllable  and 
observable,  then  the  filter  is  asymptotically  globally 
stable.  For  time  invariant  system  models,  we  would  only 
need  detectable  and  stabilizable  to  ensure  stability  of  the 
filter  [10] . 

When  Equation  (6-14)  is  applied  to  the  target  model,  it 
becomes : 


S  PI 


(6-16) 
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Ol  £  I 

J  »  i  -  N  ♦  1 


A2  Qd  AfiQd  ACQd 

AfiQd  B2  Qd  BCQd 

ACQd  BCQd  C2  Qd 


where 

A  -  t.  [At  -  B] 

B  *  tb  [1  -  C] 

C  »  exp(-At/Ti) 

Qd  *  the  discrete  target  dynamic  driving  noise  strength 
At  *=  ti  -  tj  (Note:  ti -tj  is  not  the  sample  period) 


Because  Qd  does  not  equal  zero,  and  the  matrix  in  Equation 
(6-16)  is  a  positive  definite  for  all  forward  running  time, 
there  will  always  be  positive  a  and  p  such  that  0<a<p<«; 
therefore,  the  target  model  is  fully  controllable. 

When  Equation  (6-15)  is  applied  to  the  target  model,  it 
becomes 


i 

al  s  I 

J  ■  i  -  N  ♦  1 


1/Rt  -At/Rr  C/Rt 

-At/Rr  At2  /Rt  -AtC/Rr 
C/Rt  -AtC/Rt  C2  /Rt 


*  PI  (6-17) 


where 

C  ■  exp(-At/T«) 

Rt  ■  target  measurement  variance 

At  *  tj  -  ti  (Note:  backwards  time) 


Because  R  does  not  equal  zero  or  infinity,  and  the  matrix  in 
Equation  (6-17)  is  positive  definite  for  all  backwards 


i 


i 
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running  time,  there  will  always  be  positive  a  and  $  such 
that  0<a<$<~;  therefore,  the  target  model  is  fully 
observable . 


Since  the  target  model  upon  which  the  Kalman  filter  is 
based  is  fully  controllable  and  fully  observable,  the  filter 
is  uniformly  asymptotically  stable  in  the  global  sense. 
Uniformly  stable  implies  that  the  system  in  stable  for  all 
time,  regardless  of  what  absolute  time  it  is.  The  same 
observations  can  be  made  by  referring  to  the  target  shaping 
filter  in  Figure  5-3.  It  is  obvious  that  model  is 
controllable  from  the  entry  point  of  the  white  noise,  and 
observable  from  the  point  of  extraction  of  the  target 
position  state. 

When  the  same  zero-input  stability  criteria  is  applied 
to  the  homogeneous  part  of  a  Kalman  filter  that  approximates 
the  Snyder-Fishman  filter,  the  Poisson  distributed  sample 
rate  does  not  violate  any  condition  which  define 
controllability  and  observability.  Note  that  the  Kalman  and 
Snyder-Fishman  filter  structures  are  almost  identical  (see 
Section  2.2) . 

When  stochastic  controllability  Equation  (6-14)  is 
applied  to  the  particle  beam  model,  the  resulting  expression 
is 


i 

aS  I  Qd  exp[-2  (tt -tj  ) /t*  ]  s  6  (6-18) 

J  ■  I  —  M  ♦  1 
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As  long  as  Q<i  does  not  equal  zero,  and  none  of  the  sample 
periods  are  allowed  to  approach  ••  (theoretically  possible 
with  the  Poisson  process) ,  then  there  shall  be  positive 
numbers  a  and  £  where  0<a<3<",  and  the  model  used  for  the 
basis  of  the  filter  design  will  be  controllable. 

Stochastic  observability  is  also  met.  This  can  be 
shown  by  applying  Equation  (6-15)  to  the  stochastic  particle 
beam  model  to  attain 

i 

ail  exp[2(ti-tj  )/tb]/R  i  3  (6-19) 

J  *  1  -  N  ♦  1 

But,  a  fundamental  observability  problem  could  arise  for  the 
longer  sample  periods,  because  as  the  length  of  the  sample 
period  grows,  the  upper  bound  must  grow  at  an  exponential 
rate.  Regardless,  the  Snyder-Fishman  filter  equations  are 
technically  asymptotically  stable  assuming  we  define  the 
range  of  the  Poisson  process  intervals  between  measurements 
as  not  to  include  infinity. 

6.4  Stability  Insights  on  the  Meer  Filter  Structure 

This  section  will  re-evaluate  the  method  in  which  the 
Meer  filter  gain  is  calculated  for  a  filter  depth  of  one, 
with  particular  emphasis  on  gaining  an  understanding  of  the 
severe  robustness  problems  associated  with  mismodeling  of 
the  beam  time  constant,  t« .  As  it  was  shown  in  Section  2.5, 
the  Meer  filter  gain,  Kn(ti),  is  related  to  the  Snyder- 
Fishman  filter  gain,  Ksr  (ti ) ,  by  the  following  expression 
for  a  depth  of  one: 


>\ 


a 


Km  (tt  )  -  pi  (ti  )Kar  (ti  ) 


(6-20) 


where  pi  (ti  )  is  the  probability  that  the  event  at  ti  was 


signal  induced  (see  Equation  2-27).  As  the  filter  residual, 


r(tt),  increases,  pi (ti )  decreases.  If  the  mismodeling  of 


tb  is  severe  enough,  the  residual  will  be  large  enough  so 


that  pi (ti )  will  approach  zero  and  drive  the  Meer  filter 


gain  to  zero  regardless  of  the  Snyder-Fishman  filter  gain. 


This  suggests  the  possibility  that  the  Meer  filter  structure 


can  severely  inhibit  the  robustness  of  the  controller,  and 


that  the  apparent  instability  might  be  in  part  caused  by  or 


at  least  aggravated  by  a  robustness  problem. 


The  concept  is  that  the  combined  effects  of  long 


sample  periods  and  mismodeling  tb  will  cause  large 


residuals.  Ideally,  a  large  residual  indicates  that  the 


event  is  most  likely  noise-induced,  and  the  event  receives  a 


probabilistic  weight  as  such.  But,  if  the  true  center  of 


the  beam  is  offset  from  the  filter's  estimate  of  the  center, 


then  the  Meer  filter  will  begin  to  weight  a  large  number  of 


signal-induced  events  as  being  noise-induced  (see  Figure 


6-4) .  It  appears  that  the  filter  could  counteract  this 


effect.  As  more  events  are  labeled  noise-induced,  the  time 


between  signal  events  will  appear  to  increase,  and  the 


filter  covariance  will  increase.  This  will  cause  the 


elemental  Snyder-Fishman  filter  to  rely  more  on  the  incoming 


.xv  r.  tt.  -r.  <- t. 


Figure  6-4.  True  Beam  Offset  from  Filter  Beam: 

Rt  =Rf  =0 . 5  cm,  r(ti)=1.0  cm 

measurements.  This  does  happen,  but  it  is  to  no  avail 
because  the  elemental  Snyder-Fishman  filter  gain  is 
multiplied  with  the  probabilistic  weight  that  it  was  a 
signal-induced  event,  pi (tt ) .  This  product  forms  the  Meer 
filter  gain.  As  the  filter  beam  center  diverges  from  the 
true  beam  center,  pi (ti )  approaches  zero  and  forces  the  Meer 
filter  gain  to  zero.  Meanwhile,  no  new  measurements  are 
incorporated,  and  the  filter  never  recovers. 

This  led  to  the  hypothesis  that  some  of  the  robustness 
can  be  recovered  by  computing  the  filter  R  as  a  function  of 
some  time-weighted  average  of  the  previous  N  residuals.  By 
increasing  the  filter  R  when  the  most  recent  residuals  are 
consistently  larger  than  originally  anticipated,  the  Meer 
filter  is  told  to  look  for  the  signal-induced  events  across 
a  larger  section  of  the  detector  (see  Figure  (6-5)).  The 


-9.0  -2.9  0.0  2.9  9.0 

Detector  Array 

Figure  6-5.  True  Bean;  Offset  from  Filter  Beam: 

Rt=0.5  cm,  Rf**1.0  cm,  r(ti)=  1.0  cm 

advantage  of  this  concept  is  that  most  of  the  signal-events 
can  be  recovered,  and  the  filter  should  have  much  better 
robustness  with  respect  to  the  mismodeling  of  the  beam  time 
constant.  The  robustness  analysis  found  that  the  controller 
is  rather  insensitive  to  increases  of  the  filter  R  (see 
Section  5.4.3),  and  although  the  larger  filter  R  will  cause 
filter  to  weight  more  noise-induced  events  as  signal- 
induced  events,  this  should  have  no  significant  effect  on 
the  RMS  tracking  error.  This  conclusion  is  drawn  from  the 
results  found  during  the  sensitivity  and  robustness  analyses 
conducted  on  the  beam  dispersion  parameter  (see  Sections 
5.3.3  and  5.4.3) . 

This  concept  was  simulated  through  a  measurement  model 
that  received  1000  events.  The  true  beam  had  a  dispersion 
parameter  Rt  of  0.5  cm,  and  was  located  at  the  center  of  the 
detector.  The  measurement  model  was  offset  from  the  true 
beam  center  by  a  prespecified  residual.  The  simulation 
included  noise-induced  events  as  specified  by  the  SNR,  and 


TABLE  6-2 


The  Filter’s  Probabilistic  Weight  that  an  Event  Was  Signal 
Induced  as  a  Function  of  Rf  and  the  Filter  Residual 

(Simulated  Results:  Sample  Population  =  1000  events) 


1.  Rf  *0 . 5 

Filter  residual: 

r  = 

Zb  “  £b 

(cm) 

Probability 

0.0 

0.5 

1.0 

1.5 

o 

• 

CM 

2.5 

3.0 

Si  (ti  ) 

.973 

.965 

.938 

.851 

.655 

.454 

.255 

l.OSpi < .99 

.889 

.754 

.506 

.371 

.124 

.019 

.007 

.  99£pi  <  .95 

.081 

.197 

.291 

.185 

.300 

.137 

.021 

. 955pi < .90 

.003 

.024 

.080 

.097 

.081 

.103 

.038 

. 90£pi  < .75 

.004 

.005 

.065 

.146 

.039 

.145 

.068 

.  75£pi < . 50 

.003 

.001 

.024 

.102 

.117 

.083 

.126 

. 50£pi <.30 

.000 

.002 

.008 

.034 

.085 

.016 

.082 

.  30£pi  £  .15 

.001 

.005 

.005 

.031 

.067 

.053 

.077 

. 15ipi  < .05 

.004 

.005 

.001 

.017 

.081 

.111 

.053 

. 05£pi £ . 00 

.015 

.016 

.020 

.017 

.106 

.333 

.528 

2.  Rf  =1 . 0 

Filter  residual: 

r  = 

N 

■ 

1 

(cm) 

Probability 

0.0 

0.5 

1.0 

1.5 

2.0 

2.5 

3.0 

Pi  (ti  ) 

.986 

.986 

.975 

.971 

.935 

.844 

.667 

l.OSpi <.99 

.970 

.958 

.849 

.558 

.434 

.201 

.047 

. 99Spi < .95 

.011 

.021 

.122 

.374 

.267 

.295 

.265 

. 95Spi < . 90 

.002 

.003 

.005 

.036 

.135 

.092 

.125 

.  90£pi  <  .75 

.004 

.003 

.002 

.013 

.116 

.201 

.065 

.  75ilpi  <  .  50 

.002 

.006 

.003 

.006 

.025 

.117 

.178 

. 50Spi <.30 

.001 

.001 

.001 

.000 

.007 

.048 

.140 

.30s;pi  <  .15 

.004 

.002 

.003 

.001 

.002 

.022 

.076 

. 15£pi < . 05 

.003 

.002 

.001 

.000 

.002 

.011 

.047 

. 05Spi < .00 

.003 

.004 

.014 

.012 

.012 

.013 

.054 

3.  Rf  -2 . 0 

Filter  residual: 

r  = 

Zb  -  £b 

(cm) 

Probability 

0.0 

0.5 

1.0 

1.5 

o 

• 

CM 

2.5 

3.0 

0i  (ti  ) 

.997 

.996 

.993 

.991 

.985 

.980 

.963 

l.Ospi < .99 

.976 

.980 

.980 

.962 

.834 

.598 

.476 

.  99slpi  <  .95 

.015 

.008 

.007 

.025 

.148 

.371 

.382 

.  95Spi  <  .90 

.005 

.004 

.004 

.003 

.003 

.012 

.097 

.90Spi < .75 

.003 

.006 

.003 

.003 

.004 

.007 

.031 

.  75Spi < .50 

.001 

.000 

.001 

.003 

.003 

.003 

.000 

. 50ipi < . 30 

.000 

.001 

.001 

.001 

.000 

.002 

.002 

.  30£pi  <  .15 

.000 

.001 

.002 

.001 

.002 

.000 

.003 

.  15ipi  <  .  05 

.000 

.000 

.002 

.002 

.001 

.001 

.002 

. 05£pi < . 00 

.000 

.000 

.000 

.000 

.005 

.006 

.007 

Note:  R»=0.5,  SNR=20  (952  signal-induced  events,  48  noise- 

induced  events) 


k 
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the  simulation  was  conducted  with  the  filter  R  equal  to  0.5, 
1.0  and  2.0.  The  results  are  in  Table  6-2,  and  confirm  the 
idea  that,  as  the  offset  increased,  pi (ti )  would  decrease, 
and  that  an  increase  in  R  could  compensate  for  the  offset 
and  still  be  able  to  identify  the  most  severely  degrading 
noise-induced  events  properly.  Under  nominal  conditions  and 
assuming  a  filter  residual  equal  to  0.0  cm,  the  filter 
generated  an  average  probability  that  an  event  was  signal- 
induced,  p,  of  .973.  More  significant  is  the  fact  that 
the  Meer  filter  was  able  to  distinguish  the  most  severely 
degrading  noise-induced  events  by  giving  them  an  appropriate 
probabalistic  weighting  of  .05  or  less.  As  the  filter 
residual  increased  to  3.0  cm,  p  decreased  to  .255,  and  the 
filter  had  great  difficulty  in  distinguishing  the  signal- 
induced  events.  But,  as  the  filter  beam  dispersion  (Rf )  is 
increased,  the  Meer  filter  was  able  to  identify  more  of  the 
signal-induced  events  and  up  to  some  point,  be  able  to 
distinguish  the  most  severly  degrading  noise-induced  events. 

6 . 5  Conclusions 

The  controller  stability  evaluation  led  to  several 
conclusions.  First,  the  apparent  controller  instability  is 
really  a  combination  of  an  intermittent  instability  problem 
and  a  stability  robustness  problem.  Second,  it  appears  that 
there  are  three  sources  of  the  controller  instability.  The 
varying  sample  rate  {based  on  the  mean  signal  parameter  rate 
of  one  event  per  second)  can  drive  the  controller  unstable 
during  the  longer  sample  periods.  The  mismodeling  of  tb  can 
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lead  to  large  residuals,  especially  during  the  longer  sample 
periods.  Finally,  as  the  residuals  increases,  the  Meer 
filter  structure  will  assume  that  more  of  the  incoming 
measurements  are  noise-induced  and  the  filter  will 
inadvertently  drive  the  Meer  filter  gain  to  zero.  All  three 
of  these  problems  can  be  rectified.  First,  by  insuring  the 
design  has  a  much  higher  mean  signal  rate  parameter,  we  can 
reduce  the  chance  that  a  controller  will  be  driven  unstable 
by  a  long  sample  period.  Second,  developing  on-line 
adaptive  estimation  of  tb  can  reduce  the  large  residuals. 
Third,  by  on-line  computing  of  the  filter  R  as  a  time- 
weighted  average  of  the  recent  filter  residuals,  we  can 
inhibit  the  Meer  filter  from  misdiagnosing  the  signal- 
induced  events  as  being  noise-induced,  and  thereby,  keep  the 
Meer  filter  structure  from  driving  the  filter  gain 
inappropriately  low. 


VII.  CONCLUSION  AND  RECOMMENDATIONS 


7.1  Conclusion 

Five  Monte  Carlo  analyses  and  a  PI  controller  stability 
analysis  were  performed.  The  first  Monte  Carlo  analysis 
looked  at  seven  different  cost  weighting  matrices  in  an 
effort  to  minimize  the  steady  state  RMS  tracking  error,  and 
to  enhance  the  robustness  of  the  PI  controller.  The 
significant  finding  was  that  a  10  to  1  ratio  between  the 
cost  weighting  element  on  the  states,  Xi i ,  and  the  cost 
weighting  element  on  the  pseudo-integral  term,  Xa  s  ,  resulted 
in  the  lowest  RMS  error.  This  prevents  Xa a  from  inducing 
any  degrading  oscillations  into  the  controller,  while 
allowing  the  controller  to  minimize  previous  tracking  errors 
and  compensate  for  non-linear  disturbances.  The  results  of 
next  Monte  Carlo  analysis  demonstrated  that  the  Meer  filter 
structure  can  be  reduced  to  a  depth  of  one  without  any 
apparent  degradation  of  the  controller’s  performance.  This 
reduction  in  filter  depth  should  lead  to  significant  savings 
in  on-line  computer  processing.  The  third  and  fourth  Monte 
Carlo  analyses  were  the  sensitivity  and  robustness  analyses. 
The  results  showed  the  PI  tracker  performance  was  most 
sensitive  to  equal  changes  in  the  true  and  filter  target 
dynamics  driving  noise  strengths  (Qt )  and  to  equal  changes 
in  the  true  and  filter  target  measurement  noise  variances 
(Rt ) .  This  means  that  a  PI  controller  with  perfectly  tuned 
Meer  and  Kalman  filters  (i.e.,  the  filters  have  access  to 


the  true  parameters) ,  was  most  sensitive  to  the  changes  in 
the  true  Qt  and  Rt .  When  the  controller  was  evaluated 
without  access  to  the  truth  model's  parameters,  it  was  found 
that  parameter  variations  in  the  beam  time  constant,  tb ,  and 
Qt  caused  the  most  significant  robustness  difficulties.  The 
last  Monte  Carlo  analysis  evaluated  the  on-line  adaptive 
estimation  of  Qt .  The  results  showed  that  the  multiple 
model  adaptive  controller  performed  well. 

Because  the  PI  controller  lacked  robustness  with 
respect  to  the  beam  time  constant,  tt ,  an  in-depth  stability 
analysis  was  performed  on  the  PI  controller.  This  was 
composed  of  a  zero-input  stability  analysis  on  both  the 
deterministic  full-state  feedback  and  stochastic  PI 
controllers,  a  controllability  and  observability  evaluation 
of  the  Kalman  and  Snyd.er-Fishman  filter  equations,  and  an 
analysis  on  the  Meer  structure.  The  apparent  stability 
problem  was  do  to  the  combined  effects  of  a  low  signal  rate 
parameter  that  frequently  lead  to  long,  destablizing  sample 
periods,  the  mismodeling  of  Tb  that  can  lead  to  large 
signal-induced  residuals,  and  an  inherent  design  problem 
that  would  allow  the  Meer  filter  to  evaluate  the  large 
signal-induced  residuals  as  being  noise-induced.  All  three 
problems  can  be  corrected  by  using  a  higher  signal  rate 
parameter,  by  developing  on-line  adaptive  estimation  of  tb , 
and  by  calculating  the  filter  target  measurement  noise 
variance,  Rt ,  as  a  function  of  the  time-weighted  average  of 


the  residual. 


7 . 2  Recommendations 

These  recommendations  are  divided  into  two  parts.  The 
first  part  will  suggest  ways  to  improve  SOFE  and  the  SOFE 
user-defined  subroutines  that  provide  for  the  Monte  Carlo 
analyses  in  this  thesis.  They  are  as  follows:  (1)  The 
large  number  of  runs  (200)  and  the  long  settling  time 
required  for  the  controller  (the  first  50  seconds  of  each 
run)  were  based  on  Zicker's  research  [27],  and  it  appears  to 
be  excessive.  It  is  recommend  that  the  number  of  runs  and 
the  length  of  the  initial  transient  response  be  reviewed. 
This  could  save  a  considerable  amount  of  computer  time. 

(2)  It  is  suggested  that  the  length  of  the  detector  array 
(L)  be  increased.  It  was  found  that  the  target  frequently 
exceeded  the  10  cm.  length  of  the  detector,  especially 
during  conditions  other  than  nominal.  Alternatives  would  be 
to  change  the  scales  so  that  10  centimeters  would  correspond 
to  a  large  physical  field  of  view,  or  change  the  plant  or 
target  model  gains  to  reduce  the  RMS  tracker  error. 

(3)  It  is  suggested  that  the  software  be  updated.  The  SOFE 
program  used  for  this  research  was  complied  in  May  1982  and 
was  written  in  FORTRAN  4.  The  current  version  of  SOFE  was 
written  in  FORTRAN  77  and  developed  using  top-down 
structured  programming,  and  the  current  version  (Version 
2.5)  is  supported  by  a  new  user's  manual.  The  problems  with 
using  an  older  version  of  the  software  is  the  increasing 
difficulty  in  obtaining  technical  support  and  the  inherent 
difficulty  of  reading  FORTRAN  4  code.  (The  only  SOFE  (i.e., 


167 


not  user-defined)  subroutine  altered  by  Meer  was  "ADVANS"; 
therefore,  this  should  be  the  only  subroutine  in  the  SOFE 
source  code  that  would  have  to  be  rewritten.)  Second,  if 
the  simplified  Meer  filter  (i.e.,  with  a  filter  depth  of 
one)  is  accepted  for  all  further  research,  much  of  thee 
existing  Meer  filter  code  can  be  deleted.  Most  of  the  Meer 
filter  code  is  overhead  used  to  support  greater  filter 
depths.  This  should  result  in  a  much  more  efficient  code. 

The  second  half  of  the  recommendations  suggests  ways  of 
improving  the  filter's  stability  robustness.  They  are  as 
follows:  (1)  Use  a  higher  signal  rate  parameter  that  will 

generate  a  higher  mean  sample  rate.  This  will  insure  that 
the  controller  is  stable  more  than  the  current  85  percent  of 
the  time  for  nominal  conditions.  (2)  Develop  on-line 
adaptive  estimation  of  tb ,  which  could  be  accomplished 
through  multiple  model  adaptive  methods.  If  the  Meer  filter 
is  kept  at  a  filter  depth  of  one,  this  should  be  a  simple 
task.  The  goal  is  to  reduce  the  large  residuals  that  occur 
during  the  mismodeling  of  tb .  (3)  Develop  a  variable  filter 

R  which  is  always  larger  than  the  true  R  and  is  evaluated  as 
a  function  of  the  time-weighted  average  of  the  the  most 
resent  residuals.  This  should  be  done  with  the  intent  of 
keeping  the  mean  conditional  probability  that  an  observation 
was  signal-induced,  £i (tt ) ,  near  the  true  (based  on  the 
truth  model  simulation  of  signal  and  noise  rates)  value  of 
£i  .  At  the  same  time,  the  filter  R  must  be  kept  small 
enough  to  be  able  to  identify  the  severely  degrading  noise 


events.  The  goal  is  to  maintain  stability  robustness  by 
preventing  the  Meer  filter  from  driving  pt  (tt )  to  extremely 
low  values,  which  would  effectively  "shut  down”  the  filter 
measurement  model.  (4)  Consider  using  implicit  model¬ 
following  (see  Section  3.4  and  Appendix  A)  as  a  means  to 
enhance  the  controller's  robustness  and  to  obtain  good 
response  characteristics  [16] .  Once  the  PI  controller 
demonstrates  good  robustness  for  all  meaningful  parameter 
variations,  the  effects  of  unmodeled  beam  dynamics  and 
timing  delays  can  be  studied.  (5)  Consider  reducing  the  RMS 
tracker  error  by  increasing  the  gain  of  the  particle  beam 
plant  (see  Section  5.3.1  and  specifically  Equation  (5-2)). 
This  can  be  done  be  increasing  the  control  input  gain 
matrix,  B(t)  (for  this  application,  B  is  a  scalar) .  The 
effects  of  different  control  input  gains  can  be  analyzed 
with  frequency  domain  techniques  (like  using  Bode  plots) . 
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Alternate  Cost  Weighting  Methods 


This  appendix  will  derive  the  implicit  model-following 
method  [16,19],  and  discuss  the  LQG/LTR  (linear  quadratic 
Gaussian/loop  transfer  recover)  Dual  Method  of  Kwakernaak 
and  Sivan  [10] .  Both  of  these  methods  are  designed  to 
enhance  the  controller's  stability  robustness. 

To  be  able  to  derive  the  implicit  model-following  cost 
function,  we  must  start  with  the  general  form  of  the  cost 
function: 


tN  ♦  1 

Jc  ■  EB4xT  (taw  )Xrx(tu*i  )  +  K  [xT  ( t )  W*  x  ( t )  X  ( t )  + 

,  to 


+  UT  ( t )  Wu  a  (t)u(t)  +  2xT  (t)Wxu  (t)u(t)]dt  (A-l) 


Because  we  are  designing  for  steady  state,  the  cost  function 
can  be  re-defined  as: 


••  fti  ♦  i 

Jc  “  a  l  [xT  (t)WxxX(t)  +  UT  (t)WauU(t) 

1  *°  Jtl 


+  2xT  (t)WxaU(t)  ]dt 


(A-2) 


Now,  for  all  ta[ti,ti*i),  x(t)  is  described  by: 


x(t) 


$  (t ,  ti  )x(ti  ) 
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't 

${t  , T)B(T)u(T)dT 
tl 
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“t 

0(t  , T)G(T)u(T)dP (t) 
tl 


(A-3) 


and  vi(t)  =  u(ti  )  ,  and  u(ti  )  is  constant  over  the  entire 
sample  period;  therefore  we  can  define 


B*  (t,  ti  ) 


‘t 

£>  ( t ,  T  )  B  (  T  )  dT 
tl 


(A-4) 


Therefore,  the  deterministic  form  of  Equation  (A-3)  is 

X  ( t )  *  $(t,  ti  )x(ti  )  +  BMt,  ti  )u(ti  )  (A- 5) 

and  the  last  integral  term  in  Equation  (A-3)  has  been 
dropped  because  the  only  effect  it  will  have  on  the  cost 
function  is  in  the  form  of  a  loss  function,  defined  as 
Lr (ti ) ,  and  because  it  cannot  be  affected  by  the  control 
input  [14] .  Said  another  way,  the  first  part  of  an  assumed 
certainty  equivalence  design  is  the  synthesis  of  the  LQ 
deterministic  optimal  controller,  so  the  stochastic  driving 
term  can  be  neglected.  By  substituting  Equation  (A-5)  into 
Equation  (A-2) ,  we  get 
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and  since  x(ti )  and  u(ti )  are  constants,  this  becomes 
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By  substituting  Equation  (A-6)  into 


Jc  =  H  I  [ x T  ( 1 1  )  X  ( 1 1  )x{ti  )  +  uT  (ti  )U(ti  )u(ti  ) 
i=0 

+  2xT  (ti  )  S  (ti  )u(ti  )  ]dt  ( A-7 ) 
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it  yields 
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(A-10) 


By  definition,  X(ti )  and  Wxx(t)  are  symmetrical  positive 
semi-definite,  U(ti)  and  W0B(t)  are  symmetrical  positive 
definite,  and  S(ti )  and  Wx a  (t)  are  defined  such  that  the 
expression 

(&(ti  )  -  S  (ti  )Ji -*  (ti  )S(ti  )T  ] 
is  positive  semi-definite. 

Instead  of  using  the  quadratic  cost  term 
XT  (t)  Wx  x  ( t ) 2SL ( t )  in  the  integral  in  the  integral  on  Equation 
(A-2)  that  is  used  to  drive  x(t)  to  "zero",  we  can  define  a 
model  of  desired  characteristics  as: 
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and  then  put  a  quadratic  penalty  on  the  deviations  of  the 
actual  system  outputs  (y  =  C3&)  from  those  characteristics: 


{£  -  F n y } T Wy y  (£  -  Fm yl 

»  f£*  -  Eft y I  *  W jr  y  (Cx  -  £„yl 

*  [ [CF-£mC]x  +  CBul T  Wv v  f TCF-FmCIx  +  CBu}  (A-12) 

The  jctWxbU  quadratic  cost  term  in  Equation  (A-2)  is  replaced 
with 

\±  -  FMylT  WyuU  *  {  [CF-FmC]X  +  CBu  I  T  Wy  a  U  (A~13) 

Therefore,  the  implicit  model-following  cost  function  is: 

••  fti*t 
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i »oj  t, 

+  HT  Wb.H  +  f  [££“Fm£]  S  +  Sfiuf  TWy  ■Hldt  (A-14) 

Using  the  same  derivation  as  above,  this  can  be  written  as 

•• 
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where 
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If  we  were  to  use  the  implicit  model-following  in  the 
PI  tracker  as  described  in  Section  3.3.2,  we  would  use  the 
same  augmented  states  which  include  a  pseudo-integral  term, 
C(ti )  would  be  defined  by  Equation  (3-56) , 
fi(ti)  =  [l  0  0  0  0]T  ,  and  ®(t,ti)  would  be 


0(t,ti  ) 


A  0  0  0  0 
0  1  at  D  0 
0  0  1  C  0 
0  0  0  B  0 
1-1001 


where 

A  *  exp(-At/T» ) 
B  *  exp(-At/TT) 
C  *  Tt  [1  -  B] 

D  *  Tt  [At  -  C] 
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(A-19) 


The  cost  weighting  matrices  Wxx(t),  Wxn  (t)  and  W«n  of 


the  general  cost  function,  and  Wyy  (t)  and  Wyu  (t)  of  the 
implicit  model-following  method  must  be  selected  as  to 
represent  the  physical  objective  we  are  trying  to  accomplish 
with  the  individual  states.  For  example  (see  Equation 
(3-62) ) , 


gxx  (t)  -  Wyy  (t) 


a  -a  0  0  0 

-a  a  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0000  a/10 


(A-20) 


will  drive  the  error  between  the  target  and  the  beam  states 
to  zero,  while  applying  a  10  to  1  ratio  to  the  pseudo¬ 
integral  term.  The  10  to  1  ratio  is  suggested  because  it 
delivered  the  best  performance  (see  Section  5.1).  The  cost 
weighting  matrix,  Wyy(t),  must  be  symmetrical,  positive 
semi-definite,  and  the  WyU (t)  term  must  be  such  that 

[£i  ( ti  )  -§.1  ( ti  )  U-  1  ( tt  )  Si  t  ( ti  )  ) 

is  a  positive  semi-definite  matrix. 

As  the  result  of  the  implicit  model-following 
technique,  Equation  (A-ll)  is  embedded  into  the  cost 
function,  and  the  desired  model  does  not  appear  explicitly 
in  the  final  controller  structure.  That  is,  the  additional 


states  used  to  achieve  the  desired  characteristics  of  the 
controller  are  not  augmented  to  the  original  system  model 
states.  Implicit  model-following  gives  us  the  flexibility 
to  select  the  model  of  the  desired  closed  loop  system 
characteristics  by  allowing  us  to  choose  Eu.  to  have  a 
desired  eigenstructure.  In  other  words,  we  can  select  £>  in 
such  a  manner  as  to  place  the  poles  and  orient  the 
eigenvectors  to  develop  a  controller  with  good  loop  shape 
(desirable  feedback  control  characteristics)  and  enhance 
robustness  [2] . 

Another  method  that  can  be  used  to  enhance  robustness 
is  LQG/LTR  Dual  Method  of  Kwakernaak  and  Sivan  [10] .  This 
is  a  specific  method  of  choosing  that  enhances  the 
system  stability  robustness  at  the  plant  output.  This  is 
done  by  treating  the  Kalman  filter  (or  Snyder-Fishman  filter 
as  it  would  be  the  case  for  this  application)  as  a  dynamic 
compensator  and  using  it  to  establish  good  loop  shapes  in 
order  to  achieve  good  command  following,  good  disturbance 
rejection  and  low  sensitivity  to  plant  variations  from  the 
nominal  conditions.  Then  Wxx  is  chosen  iteratively  as  some 
initial  ]frx  value  plus  the  q8£TV£.,  where  V  is  a  positive 
definite  weighting  matrix  and  q8  is  a  scaler  that  is 
increased  to  acheive  as  much  of  the  desirable  robustness 
characteristics  as  possible  in  the  implementable  filter- 
controller  structure. 


Appendix  £ 


Outline  SOFB  Operations 

B.i  g&mssarft  at  ibs.  software 

This  section  provides  a  functional  outline  of 
subroutines  used  in  SOFB  [7,21]  for  this  application.  The 
lower  case  letters  are  the  internal  SOFB  subroutines,  and 
the  upper  case  letters  are  the  user-defined  subroutines. 
Because  a  shaping  filter  was  used  to  generate  the  target 
trajectory  rather  than  using  an  externally  generated 
trajectory,  the  following  subroutines  were  not  called: 
span,  TRAJO,  move,  interp,  icsevu,  or  icsicu. 


SOFB 

***  Problem  Initialization  *** 


-yptefed 

-Rflyano 

-splita- 


-nzrcio 


-qetx  (twice) 
-getpf - zroize 


I — usm- 

— LAMBN 
— valdta- 
— plcco— 


-CNTOAIN 


-pgqrt 


-zroize 


<*> 


(Calculates  xn  from  SNR) 
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s 


[ii-qgms 

[2) -Qsrm 

[3] -fEgP9K 


if  it  was  a  signal  event,  then... 
if  it  was  a  noise  event,  then... 
if  at»l  (sec),  then  calculate  u(ti ) 


:o: 


***  update  Kalman  Piter  Equations  *** 
(The  Meer  filter  is  updated  as 
part  of  the  Propagation  cycle) 


sqrt 


- g&£ 

- XSPlUS 


- out 

***  Injects  Non-Linear  Disturbances  into  Plant  *** 

- AMEND  (called  but  not  used  for  this  thesis) 

- out 


Notes : 

GETRS ,  GETRN,  HRZ  and  SNOYS  use  the  function  GAUSS  (see 
Equations  (4-12)  and  (4-16)  for  HRZ  and  SNOYS) 

GETTS  and  GETTN  use  the  function  Poisson  (see  Equation 
(4-18) ) 
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Subroutine  "out"  calls  several  different  subroutines 
and  collects  the  data  required  for  the  SOFE  output 
file,  and  by  SOFEPL  for  statistical  reduction  of  the 
tracker  and  filter  state  errors  (see  Equations  (4-1) 
through  (4-9)). 

Definitions  of  the  Subroutines: 

advans  controls  propagation  between  measurement  updates 

advano  initializes  advans 

asysp  adds  GOGT  together  for:  dP/dt»FP+PFT  +GOGT .  the 
differental  form  of  Equation  (3-4) . 

deriv  evaluates  the  derivatives  of  the  filter  and  truth 
states 


fppft 

gauss 

getpf 

getx 

kutmer 

nzrcio 

out 

psqrt 

sof ebd 

splita 

sqrs 

valdta 

xsplus 

AMEND 

AMENDO 


adds  FP+PFT  together. 

Gaussian  random  number  generator 

reads  the  initial  covariance  matrix,  £o 

reads  in  the  initial  truth  and  filter  states,  xo 

Kutta-Merson  integration  algorithm 

reads  the  non-zero  input  from  the  disk 

outputs  all  scheduled  data 

takes  the  Cholesky  upper  triangular  square  root  of 
£(tt  ) 

reads  the  PRDATA  block  data  file  from  the  disk 
PRDATA  is  the  SOFE  simulation  specification 
file. 

partitions  the  unlabeled  common  block  'A'  into 
the  states,  filter  covariances  and  other  data 

squares  the  Cholesky  square  root  to  form  £(ti ) 

validates  the  input  data 

performs  the  measurement  update  on  the  filter  and 
truth  model. 

used  to  inject  non-linear  disturbances  into  the 
beam  dynamics  (not  used  -  variable  VTRK  set 
to  0) 

initializes  AMEND  (not  used) 
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ESTIX 


ESTIXO 

FQGEN 

FQGENO 

HRZ 

HRZO 

SNOYS 

SNOYSO 

USRIN 


INITJ 

XFDOT 

XFDOTO 

XSDOT 

XSDOTO 

GETTS 

GETTN 

LAM  BN 

GETRN 

GETRS 

GETMEAS 

SFUPDAT 


used  to  control  the  Meer  filter  updates  and 
propagation  cy jles 

used  to  initialize  the  Meer  filter  propagation 
cycle 

calculates  the  non-zero  elements  in  F(t)  and  Q(t) 
initializes  £{t)  and  £(t) 

calculates  {|t  (ti  )  ,  n  (ti  )  ,  At  (ti  )  and  zt  (ti  ) 
initializes  fiT(to) 

used  to  inject  the  white  noise  into  the  target 
shaping  filter  and  the  beam  dynamics  model 

initializes  SNOYS  (not  used) 

reads  USR1DAT  block  data  file  from  disk,  does 
initial  filter  and  gain  calculations,  sets  up 
truth  and  filter  parameter,  prints  out  non- 
scheduled  data 

initializes  the  Meer  filter  states  and  covariance 
matrices 

filter  state  derivatives 

alternate  method  for  initializing  XFDOT  (not  used) 
truth  state  derivatives 

alternate  method  for  initializing  XSDOT  (not  used) 

generates  the  time  to  the  next  signal  event 

generates  the  time  to  the  next  noise  event 

generates  the  noise  rate  parameter  (see  Equation 
(4-24) ) 

generates  the  spatial  location  of  the  noise  event, 
r*  (see  Equation  (2-2)) 

generates  the  spatial  location  of  the  signal 
event,  rs  (see  Equation  (2-1)) 

generates  the  time  to  the  next  signal,  noise  or 
measurement  update 

get  the  next  event  and  updates  the  Meer  filter 
(see  Equations  (2-3) ,  (2-13)  through  (2-15) 
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COLLAPS  implements  either  the  "best  half"  or  "merge" 
algorithms  (see  Section  2.4) 

POISSON  Poisson  random  number  generator  (see  Equation 
(4-18) ) 

SFPROP  propagates  the  Meer  filter  (see  Equations  (2-4) 
through  (2-7),  (2-10)  and  (2-11)) 

FEEDBK  calculates  u(ti )  (see  Equation  (3-63)) 

CNTGAIN  calculates  Gc *  and  g  (see  Equations  (3-64)  and 
(3-65)) 

B.2  Subroutines  Modified  to  Implement  the  MMAC 

FQGEN,  FQGENO ,  HRZ ,  HRZO,  XFDOT  and  XSDOT  were 
modified  to  included  the  increase  in  dimensionality  of 
XTf(tl),  XTf(to),  XTt(ti),  £Tt(to),  Ft  (t)  ,  Qt  (t)  ,  Pt  ( ti  )  , 

Pt (to)  and  Hr  (ti ) ,  which  occurred  as  a  result  of  embedding 
the  3  elemental  controllers  into  SOFE.  The  measurement 
updates  were  done  recursively,  using  the  same  measurement 
for  all  three  elemental  filters.  The  elemental  filter 
residuals,  rn (ti ) ,  and  residual  variances.  An (ti )  (see 
Equation  (3-80)),  were  calucated  in  subroutine  HRZ,  and  the 
probability  that  an  elemental  filter  was  correct,  pu (ti ) 
(see  Equation  (3-79)),  was  calculated  in  a  new  subroutine 
called  CPROB .  Because  the  control  input  was  the 
probabalistic  sum  of  the  elemental  control  inputs  (see 
Equation  (3-82)),  subroutine  FEEDBACK  was  replaced  with  a 
new  subroutine  called  MMAC.  CNTGAIN  was  modified  to 
calculate  the  elemental  steady-state  controller  gains, 
2.ck*(&k),  and  they  were  a  function  of  the  discretized 
uncertain  parameter,  &  .  SNOYS  was  modified  to  extend  the 


white  noise  to  all  the  elemental  filter/controllers,  and 
USRIN  had  to  be  modified  to  accommodate  the  additional 
filter/controllers.  The  parameter  variations  were  simulated 
through  subroutine  ESTIX,  and  initialized  in  ESTIXO. 


At 

pole (1) 

pole (2) 

pole (3) 

zero(l) 

zero (2) 

.01 

.99731 

.81799 

-.21466 

.82532 

.78471 

.5 

. 84321± j 

.017096 

-.26176 

.78703 

.78471 

1.0 

.74763 

.79959 

-.32636 

.75103 

.78471 

2.0 

.48560 

.81479 

-.52990 

.68687 

.78471 

3.0 

.31488 

.81769 

-.85387 

.63136 

.78471 

3.5 

.25903 

.81862 

-1.05451 

.60630 

.78471 

4.605 

.17996 

.82017 

-1.55394 

.55625 

.78471 

5.0 

.16120 

.82060 

-1.74383 

.53992 

.78471 

10.0 

.06537 

.82386 

-4.19311 

.38312 

.78471 

The  PI 

controller 

goes  unstable  at  at>3 

.369. 

C.2  Pole/Zero  Locations  for 

At  poled)  pole  (2) 

T»t -19.6 

pole (3) 

zero  (1) 

zero (2) 

.01 

.99732 

.81799 

-.21465 

.82532 

.78471 

.5 

. 84267± j 

.016724 

-.26175 

.78703 

.78471 

1.0 

.74969 

.79878 

-.32632 

.75103 

.78471 

2.0 

.48742 

.81461 

-.52973 

.68687 

.78471 

3.0 

.31649 

.81752 

-.85382 

.63136 

.78471 

3.5 

.26049 

.81845 

-1.05478 

.60630 

.78471 

4.605 

.18113 

.81999 

-1.55572 

.55625 

.78471 

5.0 

.16229 

.82041 

-1.74642 

.53992 

.78471 

10.0 

.06609 

.82364 

-4.21612 

.38312 

.78471 

The  PI  controller  goes  unstable  at  At>3.369 


C.3 


At 

pole (1) 

pole (2) 

pole (3) 

zero (1) 

zero (2) 

0.0 

1.0 

.81786 

-.21466 

.82613 

.78471 

.01 

.99733 

.81799 

-.21466 

.82532 

.78471 

.1 

.97299 

.81930 

-.22234 

- 

.78471 

.25 

.93059 

.82244 

-.23606 

- 

.78471 

.50 

. 84290±j 

.01648 

-.26175 

.78704 

.78471 

.75 

. 80872± j 

.03280 

-.29157 

- 

.78471 

1.0 

.75103 

.79822 

-.32629 

.75103 

.78471 

2.0 

.48858 

.81449 

-.52963 

.68687 

.78471 

3.0 

.31712 

.81741 

-.85378 

.63136 

.78471 

3.5 

.26142 

.81834 

-1.05494 

.60630 

.78471 

4.0 

.21939 

.81911 

-1.27438 

.58280 

.78471 

4.605 

.18167 

.81988 

-1.55878 

.55625 

.78471 

5.0 

.16299 

.82030 

-1.74806 

.53992 

.78471 

10.0 

.06655 

.82351 

-4.23078 

.38312 

.78471 

The  PI 

controller 

goes  unstable  at  at>3 

.368. 

C.4  Pole/Zero  Locations  for 

At  pole(l)  pole (2) 

Tb  t  =20 . 4 

pole (3) 

zero (1) 

zero (2) 

.01 

.99735 

.81799 

-.21466 

.82532 

.78471 

.5 

.843111 j 

.01624 

-.26174 

.78704 

.78471 

1.0 

.75235 

.79766 

-.32627 

.75103 

.78471 

2.0 

.48970 

.81438 

-.52953 

.68687 

.78471 

3.0 

.31851 

.81730 

-.85375 

.63136 

.78471 

3.5 

.26231 

.81824 

-1.05510 

.60630 

.78471 

4.605 

.18260 

.81976 

-1.55793 

.55625 

.78471 

5.0 

.16367 

.82019 

-1.74964 

.53992 

.78471 

10.0 

.06699 

.82337 

-4.24494 

.38312 

.78471 

The  PI  controller  goes  unstable  at  At>3.368 


At 

pole (1) 

pole (2) 

pole (3) 

zero ( 1 ) 

zero (2) 

.01 

.99736 

.81799 

-.21466 

.82532 

.78471 

.5 

. 84343± j 

.01588 

-.26174 

.78704 

.78471 

1.0 

.75431 

.79679 

-.32624 

.75103 

.78471 

2.0 

.49130 

.81422 

-.52939 

.68687 

.78471 

3.0 

.31992 

.81714 

-.85370 

.63136 

.78471 

3.5 

.26359 

.81808 

-1.05532 

.60630 

.78471 

4.605 

.18363 

.81960 

-1.55947 

.55625 

.78471 

5.0 

.16464 

.82002 

-1.75189 

.53992 

.78471 

10.0 

.06762 

.82318 

-4.26528 

.38312 

.78471 

The  PI 

controller 

goes  unstable  at  At>3. 

368. 

C.6 


Transfer 


-I* 


eterministic  Full-State  Feedback  Controllej 


1  Derived  in  the  s-Domain 


Yc  (s)  s  +  Gc»* 

Yr  (s)  s(s  +  1/Tit  )  +  Gc  1  *  ( S  +  Gc2*/Gcl*  ) 


(C-l) 


Assuming  we  want  to  know  the  location  of  the  poles  for 
Gct*s1.2422  and  Gc 2  * =0 . 27436 ,  Equation  (C-l)  becomes 


Yc  (s)  s  +  .27436 

xe  ^ ______ _ 

Yr  (s)  (s  +  1.024365)  (s  ♦  .267835) 


(C-2) 


These  equations  are  the  continuous-time ,  s-domain  form  of 
Equations  (6-8)  and  (6-9) . 


Appendix  0 


Sample  Data  Used  to  Calculate  the  Statistics  in  Table  5-25 

In  an  effort  to  estimate  the  time  it  takes  the  MMAC  to 
adapt  to  changes  in  Qt ,  a  single  Monte  Carlo  simulation  run 
of  2000  seconds  was  performed.  The  uncertain  parameter,  Qt , 
was  varied  from  0.01  to  1.00  cm* /sec9 ,  and  the  following 
data  was  collected. 

Qt  Increasing  Qt  Decreasing 


Sample 

to 

from 

t« 

.01  to  1.0 

ts  tT 

to 

from  1 
ts 

.  0  to 
ts 

.01 

tT 

1 

4 

1 

1 

5 

7 

34 

21 

41 

2 

9 

7 

6 

16 

6 

16 

8 

22 

3 

4 

2 

1 

6 

3 

20 

8 

23 

4 

4 

2 

1 

6 

6 

17 

4 

23 

5 

3 

12 

12 

15 

3 

33 

26 

36 

6 

7 

2 

2 

9 

4 

38 

31 

42 

7 

4 

16 

16 

20 

6 

62 

52 

68 

8 

7 

1 

1 

8 

5 

32 

13 

37 

9 

4 

8 

8 

11 

3 

34 

27 

37 

10 

3 

11 

11 

15 

3 

71 

66 

74 

11 

5 

19 

18 

22 

2 

21 

10 

23 

12 

5 

22 

22 

27 

4 

36 

26 

40 

13 

5 

7 

6 

12 

5 

31 

18 

36 

14 

6 

10 

2 

16 

3 

8 

3 

11 

15 

8 

2 

2 

10 

8 

25 

15 

33 

16 

1 

30 

14 

31 

8 

24 

17 

32 

17 

5 

2 

1 

7 

10 

23 

10 

33 

18 

9 

2 

1 

11 

5 

46 

18 

51 

19 

6 

11 

11 

17 

5 

27 

17 

32 

20 

9 

2 

2 

11 

5 

*■“ 

*■ * 

to  ,  tt  , 

ts  and  tT 

are 

defined 

in  Table 

5-25. 

Times  are  in  seconds  (or  number  of  sample  periods,  at*1.0) 
See  Table  5-25  for  more  detail. 
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Abstract 

The  goal  of  this  research  was  to '‘develop:  a  realizable 
proportional -plus-integral  (PI)  feedback  tracker  to  control 
a  neutral  particle  beam.  The  design  Is  based  on  detecting 
the  photo-electron  events  that  are  emitted  from  a  laser- 
excited  particle  beam  and  the  observed  events  are  used  by  a 
Hear  filter  to  locate  the  beam's  centerline.  The  observed 
events  are  modeled  by  a  Poisson  spacer time  process  and  are 
composed  of  both  signal r  and  noise^induced  events.  The  Meer 
filter  la  a  stochastic  multiple  model  adaptive  estimator 
which  is  composed  of  a  bank  of  Snyder  ^Fishman  filters  and  is 
designed  to  distinguish  the  signal-induced  events  from  the 
noise-induced  events.  A  target  model  is  developed  from  a 
Gauss-Harkov  acceleration  process,  and  the  target  states  are 
estimated  by  a  Kalman  filter. j  The  "optimal"  PI  controller 
design  is  based  on  th£_JJjae*x^quadratlc  (LO)  controller 
synthesis  techniqu4r'"and  the  "assumed"  certainty  equivalence 
property,  and  tWm  Kalman  filter  provides  the  reference 
(target)  stated  while  the  Meer  filter  supplies  controlled 
(beam)  states.  >The  objectives  of  the  research  were  to  (1) 
select  the  ^est"*  cost  weighting  matrices  that  minimize  the 
RHS  tracker  error  and  enhance  robustness,  (2)  simplify  the 
Heer  filter  for  easier  on-line  usage,  (3)  complete  full-'v 
scale  sensitivity  and  robustness  analyses  over  all  the  , 
Kalman  and  Heer  filter  parameters,  and  (4)  develop  on-line 
adaptive  estimation  of  those  parameters  that  great1-'  affect 
stability  robustness  and  tracker  performance.  Dur^  the 
research,  an  apparent  stability  problem  was  uncovered,  -nd  a 
fifth  objective  was  to  identify  the  source  of  the 
instability,  and  to  propose  a  solution  that  would  insure 
stability  during  parameter  variations. 


